{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Regularized OT with generic solver\n",
    "\n",
    "\n",
    "Illustrates the use of the generic solver for regularized OT with\n",
    "user-designed regularization term. It uses Conditional gradient as in [6] and\n",
    "generalized Conditional Gradient as proposed in [5][7].\n",
    "\n",
    "\n",
    "[5] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, Optimal Transport for\n",
    "Domain Adaptation, in IEEE Transactions on Pattern Analysis and Machine\n",
    "Intelligence , vol.PP, no.99, pp.1-1.\n",
    "\n",
    "[6] Ferradans, S., Papadakis, N., Peyré, G., & Aujol, J. F. (2014).\n",
    "Regularized discrete optimal transport. SIAM Journal on Imaging Sciences,\n",
    "7(3), 1853-1882.\n",
    "\n",
    "[7] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). Generalized\n",
    "conditional gradient: analysis of convergence and applications.\n",
    "arXiv preprint arXiv:1510.06567.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot\n",
    "import ot.plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data\n",
    "-------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%% parameters\n",
    "\n",
    "n = 100  # nb bins\n",
    "\n",
    "# bin positions\n",
    "x = np.arange(n, dtype=np.float64)\n",
    "\n",
    "# Gaussian distributions\n",
    "a = ot.datasets.make_1D_gauss(n, m=20, s=5)  # m= mean, s= std\n",
    "b = ot.datasets.make_1D_gauss(n, m=60, s=10)\n",
    "\n",
    "# loss matrix\n",
    "M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))\n",
    "M /= M.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VdW5x/Hvm4RBJgFBQECj4nAVUDEtiENVrFqHi7dYxeoVby04D6h1bK1tr9apiperVipaB7SKVmmpw1XUakXQUBwYZBBEQTGACDIIJFn3j3eniRDIdJJ1ht/nefZzcvbZ55w3B/LLytprr2UhBEREpOnlxS5ARCRXKYBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsEiWMLNDzWxO7Dqk9hTAIgkzO8vMPjCzdWa21MzuNbP2yWO/N7M1ybbRzDZVuf98E9QWzKzXto4JIbwRQtirAe8x1MymmtlaMytJvj7fzCx53MzsFjNbkWy3VDwm9aMAFgHM7HLgFuBnwPbAAGAX4CUzax5CODeE0CaE0Aa4CXii4n4I4QfxKndmVtDA518O3AXcBnQFugDnAgcDzZPDRgAnAfsBfYETgXMa8r65TgEsOc/M2gG/Ai4KIbwQQtgUQvgYOAUoBM6ox2sebmaLzezKpDX5uZmdZGbHmdlcM/vSzK6tcvx3zewtM/sqOfZ/zax58tjryWHvJS3uU6u8/lVmthR4sGJf8pzdk/fol9zfycyWmdnh1dS6PfBr4PwQwlMhhK+Dmx5COD2EsCE5dBjwuxDC4hDCEuB3wFl1/WykkgJYBAYCLYE/V90ZQlgDPAd8v56v2zV53e7A9cAf8DA/EDgU+IWZ7ZocWwaMBDoBBwGDgPOTOg5LjtkvaXE/UeX1O+It9RGb1f4RcBXwqJm1Ah4EHgohvFZNnQcBLYAJNXw/+wLvVbn/XrJP6kkBLOKhtzyEUFrNY58nj9fHJuDGEMIm4E/J69yVtDBnArPwP+cJIUwLIUwJIZQmre/7gO/V8PrlwC9DCBtCCOs3fzCE8AdgPjAV6AZct5XX2eL7N7PJSWt8vZlV/AJoA6yq8rxVQBv1A9efAlgElgOdttKP2i15vD5WhBDKkq8rAvKLKo+vx0MNM9vTzCYmJ/9W4/3MNQX/shDCNzUc8wegNzC6SlfCFnWy2fcfQhgYQmifPFaRE2uAdlWe1w5YEzSjV70pgEXgLWAD8MOqO82sDfADYFIT1HAv8CGwRwihHXAtUFPLcpvBl9Q/ChgL3GBmHbdyaMX3P7iG95tJ0mJP7Jfsk3pSAEvOCyGswk/CjTazY82smZkVAk8Ci4FHmqCMtsBqYI2Z7Q2ct9njXwC71fE17wKKQwg/Bf4G/L66g0IIX+Hf/z1mdrKZtTWzPDPbH2hd5dCHgcvMrLuZ7QRcDvyxjjVJFQ0auiKSLUIIt5rZCuB2YHc8DJ8FTt/Gn+6pdAUwBrgSmA48ARxZ5fEbgIfMbDv8hFvJtl7MzAYDxwJ9kl2XAe+a2ekhhHGbH598/0uS938YWAsswE/kTU4Ouw//JfBBcv/+ZJ/Uk6n7RkQkDnVBiIhEogAWEYlEASwiEokCWEQkEo2CEAA6deoUCgsLY5chklGmTZu2PITQub7PVwALAIWFhRQXF8cuQySjmNmihjxfXRAiIpEogEVyyapVMHkyLFoEugYgOgWwSLYLAcaPh969oX17OPhgKCyEHXeEq6+Gr7+OXWHOUgCLZLOlS2HQIDjlFMjPhxtvhAkT4J574Igj4JZbYK+94MUXY1eak3QSTiRbLV7s4bt4Mdx9N5xzjodwhfPOg6lTYfhwOPFEePxxGDIkXr05SC1gkWy0ZAkceqi3gP/v/+D8878dvhX694c33oDvfMdbyePHN32tOUwBLJJtNm3yMF2+HCZN8j7fbdl+e++COOggGDYMZsxomjpFASySda680kc6jB0LRUW1e06bNvDUUx7GQ4bA6tWNW6MACmCR7DJxIowaBRdf7K3guujaFZ54Aj76CC64oHHqk29RAItkizVrvK+3d2+47bb6vcZhh8G118Kjj8JLL6W2PtmCAlgkW9xwA3z6Kdx3HzRvXv/XufZa2GMPD/NvalrzUxpCASySDd5/37sehg+HgQMb9lotW/o44fnz4be/TU19Ui0FsEg2uPJKaNcObr45Na931FFw6qnelfHZZ6l5TdmCAlgk0736qg8ju/Za6Li1lefr4aaboLQUfv3r1L2mfIsCWCSTheDzOfToARdemNrX3m03v3ru/vth7tzUvrYACmCRzPbMM/D22/CrX3nfbar9/Of+ur/4RepfWxTAIhkrBPjv//YRC2ee2Tjv0aULXHqpX6I8e3bjvEcOUwCLZKrnn4fp0+Gaa6CgEefVuvRS2G671J3gk39RAItkoorW7847wxlnNO57deoEI0bAuHGwYEHjvleOUQCLZKK//x3eesuHnzVr1vjvd8UVPpvarbc2/nvlEAWwSCa69VZf0eInP2ma9+ve3WdK++MfoaSkad4zByiARTLNrFne/3vhhd4321Quuww2bPCr5CQlFMAimWbUKB8adu65Tfu+e+8Nxx/vAbx+fdO+d5ZSAItkkpISePhhH3bWuXPTv//ll8OyZX5CThpMASySSX7/e+8GGDkyzvsffjjsvz/ccYeWtU8BBbBIpti4Ee69F37wA+8OiMHMw3/2bF/uSBpEASySKZ56yhfZvOiiuHWceqp3f4weHbeOLKAAFskUo0f7ZcfHHBO3jhYt/MKMv/4VFi6MW0uGUwCLZILiYpgyxYee5aXBj+1553kdGpLWIGnwLykiNRo92lcuPuus2JW47t199eT774d162JXk7EUwCLpbvlyX634zDN91Yt0ceGF8NVX8PjjsSvJWApgkXQ3dqwPPTv//NiVfNshh0CfPnD33RqSVk8KYJF0VlbmQ88OPxz23Td2Nd9m5r8Upk/3/mmpMwWwSDp77jlYtAguuCB2JdU74wzvFvnf/41dSUZSAIuks7vvhp12gsGDY1dSvTZtfJa08ePhiy9iV5NxFMAi6WrePF/t+JxzmmbO3/o6/3zYtMlHREidKIBF0tW99/pSQ8OHx65k2/beGwYNgvvu82XspdYUwCLpaN06ePBBH2vbrVvsamp2wQXw6acwcWLsSjKKAlgkHT32mI+xTdeTb5s78UTo2dP7rKXWFMAi6SYED7I+fXysbSYoKPAJ4l9+GT78MHY1GUMBLJJu3nwT3n3XrzQzi11N7f30p9C8uYak1YECWCTdjB4N7dvD6afHrqRudtwRhg6Fhx6C1atjV5MRFMAi6WTJEnj6aTj7bGjdOnY1dXfRRbBmja+eLDVSAIukk9//HsrL02/eh9oqKoL+/b0borw8djVpTwEski6++cbH0p5wAuy2W+xq6u+ii/wikhdeiF1J2lMAi6SLceN8xeFLL41dScP86Ec+dvnOO2NXkvYUwCLpIAQPrL594YgjYlfTMM2b+wiOl1+GDz6IXU1aUwCLpIOXX4aZM33F4UwaerY155wD220Hd90Vu5K0pgAWSQd33gldusBpp8WuJDV22MFX8Hj0USgpiV1N2lIAi8T2wQfw/PN+2XGLFrGrSZ2RI2HjRi1fvw0KYJHYbr3Vx/xmyrwPtbXXXj6P8d13+9hg2YICWCSmRYt8Ucvhw6Fjx9jVpN5VV8HKlfCHP8SuJC0pgEViuuMOP+l22WWxK2kcAwbA977n3+fGjbGrSTsKYJFYSkq8ZXj66T6VY7a66ipYvBgeeSR2JWlHASwSy223+XLz11wTu5LGdeyxcOCBcOONvnSR/IsCWCSGkhK45x748Y/9ZFU2M4MbboCFC9UK3owCWCSG22/3uR9+/vPYlTSN449XK7gaCmCRprZ0qQ/NOu207G/9VqhoBS9YoKkqq1AAizS1X/3KRwTccEPsSprW8cf7qIgbbvBFR0UBLNKk5s71kQ/nnAO9esWupmmZ+UUnn32mOSISCmCRpnTttT5JzfXXx64kjkMP9RWUb74ZVqyIXU10CmCRpvL6677c0BVX+Pppueq3v/VLk3/5y9iVRKcAFmkKpaU+R+4uu8DPfha7mrj23deXXLr3Xl/9OYcpgEWawt13+6xnd94JrVrFria+3/zGp6y84IKcXjtOASzS2D77zPt8jz4aTjopdjXpoX17uOUWmDzZl7HPUQpgkcYUgo942LjRVwrOhtUuUmXYMD8pN3IkLFkSu5ooFMAijenRR2HiRLjpJthjj9jVpJe8PHjgAf/lNGKE/7LKMQpgkcayeDFcfDEcfLDfypZ69fJREc8952GcYxTAIo1h0yYYOtRHPzz4IOTnx64ofV10ka8EfdFFvjBpDlEAizSG666DN9/0q97U9bBteXkwbhy0awc/+lFOLV+kABZJtfHjfa7fc8/1VrDUrFs3eOwx+PBD+MlPcmZomgJYJJUmT4b//E8YONDH/ErtHXmkzxUxfrz/BZEDCmIXIJI15szxVYB79oQJE6Bly9gVZZ7LL4ePPvK5InbeGc47L3ZFjUoBLJIKH37oLTgzP6PfqVPsijKTGYwe7SNIzj/fT16OGBG7qkajLgiRhpoxw8/il5fDa6/ppFtDFRTAU0/5/MHnnOMXsGQpBbBIQzz/vPf3msGrr8I++8SuKDu0aOEzxw0e7MPTRo6EsrLYVaWcAlikPsrK/Oq2E07wiwnefhv+7d9iV5VdKkL40kth1Cg47jhfzimLKIBF6mrBAu9yuO46H7f6+uvQo0fsqrJTfr6PJhkzxj/nPn3g2WdjV5UyCmCR2lq71mc122cfn8f24Yfh8cehTZvYlWW/4cNh2jT/Rfcf/+H9w3Pnxq6qwRTAIjVZtcqHRe26q89j+8MfwqxZPt5Xs5s1nX32galT4fbb4Y03/P6ZZ/q/RYZSAItUp7QUXnkFzjrLr9K65ho48EC/vPixx9TlEEvz5j5WeO5cuOQS7yPed1847DBf7n7lytgV1omFHJwCTrZUVFQUiouLY5cRT1mZX0jx5ps+muHFF+HLL7174cc/9uFQ/frFrlI2t2yZz6I2dizMm+d9xocdBoMG+W2/ftC6daO9vZlNCyEU1fv5CmCBHAjgjRu9dbRsGXzxhQ/0X7QI5s/34J0xA9at82O7doWjjvK+xmOOadQfYEmREOCdd+CZZ+Bvf/Pln8C7iPbc00eo7LEHFBb6lYrduvnCqB07+r9vPbuSFMCSEkWdO4fiwYMb/4229v+tYn/Vx0PYcisv99uyssqttNS3TZtgwwbf1q/3bc0a+Ppr/7o6PXrAXntB797eWurf339g1beb2Vas8Hk5pk/3E6Zz5vgv240btzy2oADatvWtVSvYbjsfAteihXd5NGvmx+TnV2477QR33KEAltQoat48FDfVUulbC7eK/VUfN/v2lpfnt1V/GPLz/YekWbPKH5zttvOtbVvvRmjf3rfOnaFLFw/eHj38WMkN5eX+188nn/htSYn/VbRypf+S/vpr/yto/frKX+QbN/ov9tLSyl/45eW+uvWLLzY4gDUXhLi+fSGbuyBE8vK866Fbt9iV/ItGQYiIRKIAFhGJRH3AAoCZfQ3MiV1HDToBy2MXUYNMqBEyo85MqHGvEELb+j5ZfcBSYU5DTiY0BTMrVo2pkQl1ZkqNDXm+uiBERCJRAIuIRKIAlgpjYhdQC6oxdTKhzqyvUSfhREQiUQtYRCQSBbCISCQK4BxnZsea2Rwzm29mV8eup4KZ9TSzV81slpnNNLNLkv0dzewlM5uX3HZIg1rzzWy6mU1M7u9qZlOTz/QJM2seub72ZvaUmX1oZrPN7KB0+xzNbGTy7zzDzB43s5bp8Dma2QNmVmJmM6rsq/azM/c/Sb3vm1mN85cqgHOYmeUDdwM/APYBTjOzdFnWtxS4PISwDzAAuCCp7WpgUghhD2BScj+2S4DZVe7fAtwZQugFrATOjlJVpbuAF0IIewP74bWmzedoZt2Bi4GiEEJvIB8YSnp8jn8Ejt1s39Y+ux8AeyTbCODeGl89hKAtRzfgIODFKvevAa6JXddWap0AfB+/Wq9bsq8bfgFJzLp6JD+ERwITAcOv3iqo7jOOUN/2wEKSE+5V9qfN5wh0Bz4FOuIXh00EjkmXzxEoBGbU9NkB9wGnVXfc1ja1gHNbxX/8CouTfWnFzAqBA4CpQJcQwufJQ0uBLpHKqjAKuBIoT+7vAHwVQihN7sf+THcFlgEPJt0k95tZa9LocwwhLAFuBz4BPgdWAdNIr8+xqq19dnX+eVIAS1ozszbA08ClIYTVVR8L3syINo7SzE4ASkII02LVUAsFQD/g3hDCAcBaNutuSIPPsQMwGP9lsRPQmi3/7E9LDf3sFMC5bQnQs8r9Hsm+tGBmzfDwHRdC+HOy+wsz65Y83g0oiVUfcDDw72b2MfAnvBviLqC9mVXMsxL7M10MLA4hTE3uP4UHcjp9jkcBC0MIy0IIm4A/459tOn2OVW3ts6vzz5MCOLe9A+yRnG1ujp/4+EvkmgA/owyMBWaHEO6o8tBfgGHJ18PwvuEoQgjXhBB6hBAK8c/ulRDC6cCrwMnJYbFrXAp8amZ7JbsGAbNIo88R73oYYGatkn/3ihrT5nPczNY+u78AZyajIQYAq6p0VVQvVse7tvTYgOOAucBHwHWx66lS1yH4n3bvA+8m23F4H+skYB7wMtAxdq1JvYcDE5OvdwPeBuYD44EWkWvbHyhOPstngQ7p9jkCvwI+BGYAjwAt0uFzBB7H+6U34X9NnL21zw4/AXt38rP0AT6qY5uv36BLkc3sWPxPrnzg/hDCzfV+MRGRHFPvAE7GkM7FhwYtxv+cPS2EMCt15YmIZK+GTMj+XWB+CGEBgJn9CT+TudUA7tSpUygsLGzAW0pDffWVL/7as+e390+bNm15CKFzxf3v5/1IszSJVPFS+fitLOddfw0J4OrGvPXf/CAzG4FfFcLOO+9MsVbejerKK2H06C0XQDazRXEqEsldjT4KIoQwJoRQFEIo6ty5c81PkEa1aRM0axa7ChGBhgVwWo8hlept3AjNo04NIyIVGhLAaTuGVLZu3Tpo3Tp2FSICDegDDiGUmtmFwIv4MLQHQggzU1aZNIrVq6FNm9hViAg0cFn6EMJzwHMpqkWawIoVsMMOsasQEdClyDlnyRLo2jV2FSICCuCcUlYGH38Mu+8eu5KIzHwTSQMK4BwyYwaUlkLv3rErERFoYB+wZJa//91vDzssbh1RVVx6X9EKtqQNEsq//bhIE1ALOIc8+ijss8+WlyGLSBxqAeeIKVPgnXf8MmShsqUbyvw2Lx8Aa+4/EmFTshJOeVlTVyY5RC3gHLBxI5x3HnTpAmeeGbsaEamgFnCWCwGuuw7efReefRbatYtdUZpKWrphg99agf9o2HZ+2WDYuMlvN22MUJxkK7WAs1gIPvvZ7bd7C3jw4NgViUhVagFnqVWr4OKL4eGH4YIL4H/+J3ZFmSWUln7rNq9lS7/d3i8jDOvWA1C+bl2E6iRbqAWchZ591kc7PPooXH+9n3jL07+0SNpRCziLvPUW/Pa38Ne/Qt++HsTf+U7sqrJD+Tff+BfJbX7SmV6wi4/pC6vXAFC2cmXTFycZS+2iDFdWBs88AwcfDAMHwj/+ATfd5CteKHxF0ptawBlq9mx44gkYNw7mz4fCQu/n/a//0nSTTaFs9Wr/IrnNT1Z7yeu7t98u+wqA0s+XNn1xkjEUwBlk7lwP3Sef9HkdzPyy4htvhB/+EAr0rymSUfQjm8bWrIE33oCXX4aXXoIPPvDQPeQQP7E2ZAh06xa7SgEoW7bMv0hu8wp39v1H9AOg+WJvEZfNW9D0xUnaUgCnkU2bYOpUD9xJk/zy4dJSX8Pt4INh1Cg4+WTo3j12pSKSCgrgiEpKPHArtrfegrVrvZV74IFw+eVw1FEevtttF7taqYvSjz8BID+5Zd+9APh66AAA2s37GoAwTat45TIFcBPZsMEvB54yxcN2yhRYuNAfy8+HPn1g2DAYNAgOPxw6doxarog0AQVwI1izBt5/3wN3+nS/ff99nxQHvAthwAC/PLh/f2/taqXi7FY2cw4AbZMGb/nA/QD4YuRAADq9twGAglemNX1xEo0CuIG++OLbQTt9OsybVznbYceOcMABcMklHrb9+0OPHnFrFpH0oACupfXrYdYsH4lQdVtaZZhnYaGH7emnw/77+9c9emgJMtmSTX4PgK6T/f6G4/yqmUW3HQRAt7d8hY5Wf57a9MVJk1EAb6asDBYs2DJo58+H8mTVmpYtfa6FY46pDNr99oP27ePWLiKZJacDuKRky6CdORMqJrgy8xWE+/SBoUP9tk8f6NXLT5yJpEqL594BYPfn/P7qH/toiYV/6gvA9v/nJwk6PvBW0xcnjSYnAnjdOg/WzcO2pKTymM6dPVyHD68M2n331ckxEWk8WRfAX37pJ8L++c/K27lzK0+KtWzpy7Iff3xl0Pbp48v1iKSLdo9NSW79fsn5Plqi9es+58RHz+wBQNc7Jzd9cZIyGRvAIcBnn20Ztp98UnlMz57Qr9+3uw92313dByKSHjImgMvLvRvh9dd9foTXX4fPP698fM894aCDfPWHAw7wrVOnePWKpNKO93hLd+09fr/sam8BnzDT5x++90/HA9DzN2oRZ5K0DeBNm2DatMqwffNNqJjrunt3v1pswABv4e63H7RtG7VcEZE6S6sADsEnFH/4YZ9ysWLK1T339OkWDz3Up18sLNTYWslt3W/2lu7EmzsAUPYbP8lxy0IfN3zKE5cCsOvVGjWRzmoMYDPrCTwMdAECMCaEcJeZdQSeAAqBj4FTQgj1Wo9lwQJ45BEP3gULfOTBkCFw4ok+9WLXrvV5VRGR9FabFnApcHkI4Z9m1haYZmYvAWcBk0IIN5vZ1cDVwFV1LeB3v4MrrvAW7aBBcMMN3trV8C+R2iv8hbd0r/pFfwDKb/P9Ty/20RT7Jy3i3S+f0vTFyVbVGMAhhM+Bz5Ovvzaz2UB3YDBweHLYQ8Br1DGAx4zx8B0yBO6800ctiIjkijr1AZtZIXAAMBXokoQzwFK8i6LWFi6Ec8+FoiJ47DGfdFxEUmP3n3mLeMjP/Iq6cKfvf/Gzd/3xP50LQK/L1CKOqdarIptZG+Bp4NIQwuqqj4UQAt4/XN3zRphZsZkVL6tYtgWfpOaII3z87nPP1a94EZFMVqsANrNmePiOCyH8Odn9hZl1Sx7vBpRU99wQwpgQQlEIoahzsnIsQLNm8Oyz3gIeMgSOO84XnPzmmwZ9PyJSjV4jp9Br5BSO2Wl/jtlpfyyABe8jfnrxFD667SA+SmZik6ZTYwCbmQFjgdkhhDuqPPQXYFjy9TBgQl3fvG1beOEFuOoqn5th6FAf8TBihI/7rZh9TEQkG1kI1fYcVB5gdgjwBvABUBGJ1+L9wE8COwOL8GFoX27rtYqKikJxcXG1j5WVwWuvwUMPwdNP+wQ6HTr4MLTDDvMxwP36ectZUs/MpoUQiirufz/vR9v+jyFZZeHN3vp98tRRAJz6qI+aqBhdIfBS+fiUX31Qm1EQ/wC29saDUlVIfr4PQxs0CO65ByZMgFdf9Svh/vpXP6ZVK7/cuOKCjO9+V8PVRCRz1dgCTqVttYC3ZelSv0KuYh6I997zq+bMYK+9fN6Hfv0q54DQgpZ1pxawVPXJ9T772vmn/g2A+x7xuSYqrsDLRVFawOmga1c4+WTfAL76CiZPhnfe8VEU//gHPP545fG77FIZyP36+aoVO+2ky5dFJL1kRABvrn17HzVx3HGV+5Yv9zCu2P75Tx9lUdHA79ChckrKvn39tndvTeIjUp2df53MNfFrn2sif6Tvr5iP+OPHewHQ+V71ETdERnRB1NfXX3t3xXvv+SiL99+HGTN8f4XCwm9PzN6nj0/+k2sn+9QFIXXx5U/8pN2qo9cCsMNfWgGVE8lno5ztgqivtm19FMUhh1TuCwEWLdpyeaLnn4fSUj+meXPYe+8tg1krHItIKmV1C7guNmyAOXO8lVw1mBcvrjymfXvvtti8K6Ndu3h1p4pawNIQa4f4JEBLB/ilBT1e9dZMxWKj2UAt4EbUooUHat++396/cqV3W1QN5XHjKucqBl8luWJ5+orbrl3VWhaRbVMLuB5CgE8/9dbyu+/6Nn26z2VcYccdKwN5//3hwAM9qNM1lNUCllQqPfJAAJbv14IuU72f2Ca/F7OkBlMLOE2Ywc47+3bCCZX7V63yUJ4+vTKU77jDl1cCH5/cv79vAwb4hSQdOsT5HkQkPgVwCm2/vV+ld+ihlfs2boRZs6C4GKZOhSlTfP6Lij889tzTw7h/fxg40LtA8mo9R51Ieip4ZRoAXV8BO3BfAFYP9akx23/gC+eUzZwTp7g0ogBuZM2bV3ZD/PSnvm/1ag/kKVM8lF94wZdjAl/J+Ygj4Kij/LLs3XZL324LEWkYBXAE7drBkUf6BpVD415/HSZN8m38eH+ssLByjoxjj1WXhWSeMG0mAG2nJTv22A2AsiP6AdBigc8TXrro0yavLTYFcBow86AtLIQzz/RAnjOnMoyffhrGjvWLQ44+Gk45BQYP9i4PEclcCuA0ZOYXguy9N1xwgU/V+c47HsRPPgl/+5t3bRx7rIfxkCHQsmXsqkVqp2yeDxfKn5fs6L4TAHl99wbAPl/hx1VZQSdb6XRPBsjP9xN1t90GH38Mb73lwTxtGpxxhk8+9JvfwIoVsSsVkbpQAGcYMw/jO+6ATz6Bl1/2McbXX++rSl94ofcni2SK0iWfUbrkM8rf/5Dy9z/0OQFKSynYpScFu/Qkr21b8rJ01iwFcAbLy/OTc889V7mk05gxsM8+MGqUd12ISPpSAGeJ3r3hgQdg3jw4/HAYOdLHFc+cGbsykbopW7mSspUrKV30qY+MKC+H8nLyO+1AfqcdyGvZkrwsOemhAM4yu+wCEyfCY4/5pdEDB/oCpyKSfhTAWcgMTjvNJ6Xv2tWHrk2aFLsqkfopX7uW8rVrKVu+grLlKwghEEIgr3Vr8lq3xgoKsILMHNClAM5iPXv6xR277ebLOS1ZErsiEalKAZz3yuHQAAAGyklEQVTlunSBZ57x+Y6HD6+cg0IkU4UNGwgbNvyrZVzBWrTAWrTwPwEz5Pp9BXAO6NULbrrJV/149dXY1YhIBQVwjjj3XJ8O8557YlciklqhtNS3pGX8L3n5vqUxBXCOaNkShg2DCROgyl9tIhKRAjiHHH20X2Q0JXsXrhXxEx0hQHmZbxV9wmnYN6wAziEH+UrivP123DpExGXm4Dmpl+2391ERH30UuxKRJpTGQ3/UAs4xhYU+iY+IxKcAzjGdOmnaSpF0oQDOMR06wMqVsasQEahDAJtZvplNN7OJyf1dzWyqmc03syfMrHnjlSmp0qoVrF8fuwoRgbq1gC8BZle5fwtwZwihF7ASODuVhUnjaNlSASySLmoVwGbWAzgeuD+5b8CRwFPJIQ8BJzVGgZJazZrBpk2xqxARqH0LeBRwJVCe3N8B+CqEUJrcXwx0r+6JZjbCzIrNrHhZDiyyl+4KCvxiDBGJr8YANrMTgJIQwrT6vEEIYUwIoSiEUNS5c+f6vISkUF6eLzAgIvHV5kKMg4F/N7PjgJZAO+AuoL2ZFSSt4B6AZpvNAApgkfRRYws4hHBNCKFHCKEQGAq8EkI4HXgVODk5bBgwodGqlJQxS+sLg0RySkPGAV8FXGZm8/E+4bGpKUkakwJYJH3UaS6IEMJrwGvJ1wuA76a+JGlMaTYZlEhO05VwIiKRKIBzjFrAIulDASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikdQqgM2svZk9ZWYfmtlsMzvIzDqa2UtmNi+57dDYxYqIZJPatoDvAl4IIewN7AfMBq4GJoUQ9gAmJfdFRKSWagxgM9seOAwYCxBC2BhC+AoYDDyUHPYQcFJjFSkiko1q0wLeFVgGPGhm083sfjNrDXQJIXyeHLMU6FLdk81shJkVm1nxsmXLUlO1iEgWqE0AFwD9gHtDCAcAa9msuyGEEIBQ3ZNDCGNCCEUhhKLOnTs3tF4RkaxRmwBeDCwOIUxN7j+FB/IXZtYNILktaZwSRUSyU40BHEJYCnxqZnsluwYBs4C/AMOSfcOACY1SoYhIliqo5XEXAePMrDmwAPgvPLyfNLOzgUXAKY1ToohIdqpVAIcQ3gWKqnloUGrLERHJHboSTkQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiERSqwA2s5FmNtPMZpjZ42bW0sx2NbOpZjbfzJ4ws+aNXayISDapMYDNrDtwMVAUQugN5ANDgVuAO0MIvYCVwNmNWaiISLapbRdEAbCdmRUArYDPgSOBp5LHHwJOSn15IiLZq8YADiEsAW4HPsGDdxUwDfgqhFCaHLYY6F7d881shJkVm1nxsmXLUlO1iEgWqE0XRAdgMLArsBPQGji2tm8QQhgTQigKIRR17ty53oWKiGSb2nRBHAUsDCEsCyFsAv4MHAy0T7okAHoASxqpRhGRrFSbAP4EGGBmrczMgEHALOBV4OTkmGHAhMYpUUQkO9WmD3gqfrLtn8AHyXPGAFcBl5nZfGAHYGwj1ikiknUKaj4EQgi/BH652e4FwHdTXpGISI7QlXAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwDlmwAC45JLYVYgIgIUQmu7NzJYBi5rsDaUudgkhdI5dhEguadIAFhGRSuqCEBGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKR/D+iA6AGU5wC6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% EMD\n",
    "\n",
    "G0 = ot.emd(a, b, M)\n",
    "\n",
    "pl.figure(3, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, G0, 'OT matrix G0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD with Frobenius norm regularization\n",
    "--------------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "    0|1.760578e-01|0.000000e+00\n",
      "    1|1.669467e-01|-5.457501e-02\n",
      "    2|1.665639e-01|-2.298130e-03\n",
      "    3|1.664378e-01|-7.572776e-04\n",
      "    4|1.664077e-01|-1.811855e-04\n",
      "    5|1.663912e-01|-9.936787e-05\n",
      "    6|1.663852e-01|-3.555826e-05\n",
      "    7|1.663814e-01|-2.305693e-05\n",
      "    8|1.663785e-01|-1.760450e-05\n",
      "    9|1.663767e-01|-1.078011e-05\n",
      "   10|1.663751e-01|-9.525192e-06\n",
      "   11|1.663737e-01|-8.396466e-06\n",
      "   12|1.663727e-01|-6.086938e-06\n",
      "   13|1.663720e-01|-4.042609e-06\n",
      "   14|1.663713e-01|-4.160914e-06\n",
      "   15|1.663707e-01|-3.823502e-06\n",
      "   16|1.663702e-01|-3.022440e-06\n",
      "   17|1.663697e-01|-3.181249e-06\n",
      "   18|1.663692e-01|-2.698532e-06\n",
      "   19|1.663687e-01|-3.258253e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   20|1.663682e-01|-2.741118e-06\n",
      "   21|1.663678e-01|-2.624135e-06\n",
      "   22|1.663673e-01|-2.645179e-06\n",
      "   23|1.663670e-01|-1.957237e-06\n",
      "   24|1.663666e-01|-2.261541e-06\n",
      "   25|1.663663e-01|-1.851305e-06\n",
      "   26|1.663660e-01|-1.942296e-06\n",
      "   27|1.663657e-01|-2.092896e-06\n",
      "   28|1.663653e-01|-1.924361e-06\n",
      "   29|1.663651e-01|-1.625455e-06\n",
      "   30|1.663648e-01|-1.641123e-06\n",
      "   31|1.663645e-01|-1.566666e-06\n",
      "   32|1.663643e-01|-1.338514e-06\n",
      "   33|1.663641e-01|-1.222711e-06\n",
      "   34|1.663639e-01|-1.221805e-06\n",
      "   35|1.663637e-01|-1.440781e-06\n",
      "   36|1.663634e-01|-1.520091e-06\n",
      "   37|1.663632e-01|-1.288193e-06\n",
      "   38|1.663630e-01|-1.123055e-06\n",
      "   39|1.663628e-01|-1.024487e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   40|1.663627e-01|-1.079606e-06\n",
      "   41|1.663625e-01|-1.172093e-06\n",
      "   42|1.663623e-01|-1.047880e-06\n",
      "   43|1.663621e-01|-1.010577e-06\n",
      "   44|1.663619e-01|-1.064438e-06\n",
      "   45|1.663618e-01|-9.882375e-07\n",
      "   46|1.663616e-01|-8.532647e-07\n",
      "   47|1.663615e-01|-9.930189e-07\n",
      "   48|1.663613e-01|-8.728955e-07\n",
      "   49|1.663612e-01|-9.524214e-07\n",
      "   50|1.663610e-01|-9.088418e-07\n",
      "   51|1.663609e-01|-7.639430e-07\n",
      "   52|1.663608e-01|-6.662611e-07\n",
      "   53|1.663607e-01|-7.133700e-07\n",
      "   54|1.663605e-01|-7.648141e-07\n",
      "   55|1.663604e-01|-6.557516e-07\n",
      "   56|1.663603e-01|-7.304213e-07\n",
      "   57|1.663602e-01|-6.353809e-07\n",
      "   58|1.663601e-01|-7.968279e-07\n",
      "   59|1.663600e-01|-6.367159e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   60|1.663599e-01|-5.610790e-07\n",
      "   61|1.663598e-01|-5.787466e-07\n",
      "   62|1.663596e-01|-6.937777e-07\n",
      "   63|1.663596e-01|-5.599432e-07\n",
      "   64|1.663595e-01|-5.813048e-07\n",
      "   65|1.663594e-01|-5.724600e-07\n",
      "   66|1.663593e-01|-6.081892e-07\n",
      "   67|1.663592e-01|-5.948732e-07\n",
      "   68|1.663591e-01|-4.941833e-07\n",
      "   69|1.663590e-01|-5.213739e-07\n",
      "   70|1.663589e-01|-5.127355e-07\n",
      "   71|1.663588e-01|-4.349251e-07\n",
      "   72|1.663588e-01|-5.007084e-07\n",
      "   73|1.663587e-01|-4.880265e-07\n",
      "   74|1.663586e-01|-4.931950e-07\n",
      "   75|1.663585e-01|-4.981309e-07\n",
      "   76|1.663584e-01|-3.952959e-07\n",
      "   77|1.663584e-01|-4.544857e-07\n",
      "   78|1.663583e-01|-4.237579e-07\n",
      "   79|1.663582e-01|-4.382386e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   80|1.663582e-01|-3.646051e-07\n",
      "   81|1.663581e-01|-4.197994e-07\n",
      "   82|1.663580e-01|-4.072764e-07\n",
      "   83|1.663580e-01|-3.994645e-07\n",
      "   84|1.663579e-01|-4.842721e-07\n",
      "   85|1.663578e-01|-3.276486e-07\n",
      "   86|1.663578e-01|-3.737346e-07\n",
      "   87|1.663577e-01|-4.282043e-07\n",
      "   88|1.663576e-01|-4.020937e-07\n",
      "   89|1.663576e-01|-3.431951e-07\n",
      "   90|1.663575e-01|-3.052335e-07\n",
      "   91|1.663575e-01|-3.500538e-07\n",
      "   92|1.663574e-01|-3.063176e-07\n",
      "   93|1.663573e-01|-3.576367e-07\n",
      "   94|1.663573e-01|-3.224681e-07\n",
      "   95|1.663572e-01|-3.673221e-07\n",
      "   96|1.663572e-01|-3.635561e-07\n",
      "   97|1.663571e-01|-3.527236e-07\n",
      "   98|1.663571e-01|-2.788548e-07\n",
      "   99|1.663570e-01|-2.727141e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  100|1.663570e-01|-3.127278e-07\n",
      "  101|1.663569e-01|-2.637504e-07\n",
      "  102|1.663569e-01|-2.922750e-07\n",
      "  103|1.663568e-01|-3.076454e-07\n",
      "  104|1.663568e-01|-2.911509e-07\n",
      "  105|1.663567e-01|-2.403398e-07\n",
      "  106|1.663567e-01|-2.439790e-07\n",
      "  107|1.663567e-01|-2.634542e-07\n",
      "  108|1.663566e-01|-2.452203e-07\n",
      "  109|1.663566e-01|-2.852991e-07\n",
      "  110|1.663565e-01|-2.165490e-07\n",
      "  111|1.663565e-01|-2.450250e-07\n",
      "  112|1.663564e-01|-2.685294e-07\n",
      "  113|1.663564e-01|-2.821800e-07\n",
      "  114|1.663564e-01|-2.237390e-07\n",
      "  115|1.663563e-01|-1.992842e-07\n",
      "  116|1.663563e-01|-2.166739e-07\n",
      "  117|1.663563e-01|-2.086064e-07\n",
      "  118|1.663562e-01|-2.435945e-07\n",
      "  119|1.663562e-01|-2.292497e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  120|1.663561e-01|-2.366209e-07\n",
      "  121|1.663561e-01|-2.138746e-07\n",
      "  122|1.663561e-01|-2.009637e-07\n",
      "  123|1.663560e-01|-2.386258e-07\n",
      "  124|1.663560e-01|-1.927442e-07\n",
      "  125|1.663560e-01|-2.081681e-07\n",
      "  126|1.663559e-01|-1.759123e-07\n",
      "  127|1.663559e-01|-1.890771e-07\n",
      "  128|1.663559e-01|-1.971315e-07\n",
      "  129|1.663558e-01|-2.101983e-07\n",
      "  130|1.663558e-01|-2.035645e-07\n",
      "  131|1.663558e-01|-1.984492e-07\n",
      "  132|1.663557e-01|-1.849064e-07\n",
      "  133|1.663557e-01|-1.795703e-07\n",
      "  134|1.663557e-01|-1.624087e-07\n",
      "  135|1.663557e-01|-1.689557e-07\n",
      "  136|1.663556e-01|-1.644308e-07\n",
      "  137|1.663556e-01|-1.618007e-07\n",
      "  138|1.663556e-01|-1.483013e-07\n",
      "  139|1.663555e-01|-1.708771e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  140|1.663555e-01|-2.013847e-07\n",
      "  141|1.663555e-01|-1.721217e-07\n",
      "  142|1.663554e-01|-2.027911e-07\n",
      "  143|1.663554e-01|-1.764565e-07\n",
      "  144|1.663554e-01|-1.677151e-07\n",
      "  145|1.663554e-01|-1.351982e-07\n",
      "  146|1.663553e-01|-1.423360e-07\n",
      "  147|1.663553e-01|-1.541112e-07\n",
      "  148|1.663553e-01|-1.491601e-07\n",
      "  149|1.663553e-01|-1.466407e-07\n",
      "  150|1.663552e-01|-1.801524e-07\n",
      "  151|1.663552e-01|-1.714107e-07\n",
      "  152|1.663552e-01|-1.491257e-07\n",
      "  153|1.663552e-01|-1.513799e-07\n",
      "  154|1.663551e-01|-1.354539e-07\n",
      "  155|1.663551e-01|-1.233818e-07\n",
      "  156|1.663551e-01|-1.576219e-07\n",
      "  157|1.663551e-01|-1.452791e-07\n",
      "  158|1.663550e-01|-1.262867e-07\n",
      "  159|1.663550e-01|-1.316379e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  160|1.663550e-01|-1.295447e-07\n",
      "  161|1.663550e-01|-1.283286e-07\n",
      "  162|1.663550e-01|-1.569222e-07\n",
      "  163|1.663549e-01|-1.172942e-07\n",
      "  164|1.663549e-01|-1.399809e-07\n",
      "  165|1.663549e-01|-1.229432e-07\n",
      "  166|1.663549e-01|-1.326191e-07\n",
      "  167|1.663548e-01|-1.209694e-07\n",
      "  168|1.663548e-01|-1.372136e-07\n",
      "  169|1.663548e-01|-1.338395e-07\n",
      "  170|1.663548e-01|-1.416497e-07\n",
      "  171|1.663548e-01|-1.298576e-07\n",
      "  172|1.663547e-01|-1.190590e-07\n",
      "  173|1.663547e-01|-1.167083e-07\n",
      "  174|1.663547e-01|-1.069425e-07\n",
      "  175|1.663547e-01|-1.217780e-07\n",
      "  176|1.663547e-01|-1.140754e-07\n",
      "  177|1.663546e-01|-1.160707e-07\n",
      "  178|1.663546e-01|-1.101798e-07\n",
      "  179|1.663546e-01|-1.114904e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  180|1.663546e-01|-1.064022e-07\n",
      "  181|1.663546e-01|-9.258231e-08\n",
      "  182|1.663546e-01|-1.213120e-07\n",
      "  183|1.663545e-01|-1.164296e-07\n",
      "  184|1.663545e-01|-1.188762e-07\n",
      "  185|1.663545e-01|-9.394153e-08\n",
      "  186|1.663545e-01|-1.028656e-07\n",
      "  187|1.663545e-01|-1.115348e-07\n",
      "  188|1.663544e-01|-9.768310e-08\n",
      "  189|1.663544e-01|-1.021806e-07\n",
      "  190|1.663544e-01|-1.086303e-07\n",
      "  191|1.663544e-01|-9.879008e-08\n",
      "  192|1.663544e-01|-1.050210e-07\n",
      "  193|1.663544e-01|-1.002463e-07\n",
      "  194|1.663543e-01|-1.062747e-07\n",
      "  195|1.663543e-01|-9.348538e-08\n",
      "  196|1.663543e-01|-7.992512e-08\n",
      "  197|1.663543e-01|-9.558020e-08\n",
      "  198|1.663543e-01|-9.993772e-08\n",
      "  199|1.663543e-01|-8.588499e-08\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  200|1.663543e-01|-8.737134e-08\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecVPW9//HXZwtdqUsvSxEQVBQ22EtEDcYeE0uM0UQv1xsTYxJ7mrk39iS2XHvD39VEo8YWS+wNQVlRUBREOtKLdLZ9fn98z7gLLGyfc2b2/Xw8zmN2Zs7MfIZh9z3f7/me79fcHRERkaTJibsAERGR6iigREQkkRRQIiKSSAooERFJJAWUiIgkkgJKREQSSQElIvVmZgeb2YwE1DHXzI6Iuw5pXAookRiZ2dlmNs3MNprZEjO73cw6RPfdYWbro63EzEqrXH8+DbW5mQ3a2T7u/pa7D6nn8881s01V3tN6M+tZv2olGymgRGJiZr8CrgMuBtoD+wH9gJfMrIW7n+fu7dy9HXA18EjqursfHV/lgZnlNcLTHFflPbVz9y+b6HVqzcxy0/l6smMKKJEYmNmuwB+An7n7C+5e6u5zgVOAQuAH9XjOw8xsoZldYmbLzGyxmZ1oZt82s5lmtsrMrqiy/2gze9fM1kT7/tXMWkT3vRnt9lHUsjm1yvNfamZLgPtTt0WPGRi9xsjoek8zW25mh9XxfRRGrbdzzGw+8Gp0+/Fm9klU7+tmtvs2D/2GmU03s9Vmdr+Ztarl6z0QtVyfM7MNwDfNrKWZ/cnM5pvZ0qg127rKYy6J/s2+NLNza9PalLpTQInE4wCgFfBE1RvdfT3wHHBkPZ+3e/S8vYDfAXcTwm4UcDDwWzPrH+1bDvwC6ALsD4wBfhLVcUi0z4ioZfNIlefvRGjpjdum9i+AS4H/M7M2wP3AeHd/vZ7v5VBgd+BbZjYY+BtwIVBA+Dd6JhWokTOAbwEDgcHAb+rwWt8HrgJ2Ad4Gro2eY29gEJX/npjZWOCXwBHRfYfV691JjRRQIvHoAqxw97Jq7lsc3V8fpcBV7l4K/D16npvdfZ27fwJMB0YAuHuxu09097Ko9XYnIRR2pgL4vbtvcfdN297p7ncDs4BJQA/g1zU835NRi2iNmT25zX1XuvuG6HVOBf7l7i9F7+1PQGtC0Kf81d0XuPsqQticXsNrV/WUu7/j7hXAFkL4/sLdV7n7OkIX62nRvqcA97v7J+6+EbiyDq8jdZDWvl0R+doKoIuZ5VUTUj2i++tjpbuXRz+nAmRplfs3Ae0AolbJX4AioA3h70FxDc+/3N0317DP3cDTwDh331LDvie6+8s7uG9BlZ97AvNSV9y9wswWEFo21e0/L3pMbVV9bAHh36PYzFK3GZA6NtUTmLyDx0ojUgtKJB7vEr6pf6fqjWbWDjgaeCUNNdwOfAbs5u67AlcQ/hDvzE6XP4jqvwm4F7jSzDo1oL6qr/UloVsx9ToG9AEWVdmnT5Wf+0aPqc9rrSAE+XB37xBt7aPBKhBauL138LrSiBRQIjFw968IgyRuNbOxZpZvZoXAo8BC4P+loYxdgLXAejMbCvzXNvcvBQbU8TlvBia7+7nAv4A7Glxl8ChwjJmNMbN84FeEgJ9QZZ/zzax3FIq/Bh6p5nlqFHXz3Q3caGZdAcysl5l9q0otPzKz3aNjbb+t31uSmiigRGLi7tcTWi1/IgTFJEJ30ZhadI01hosIgwPWEf4gb/sH/UpgfHR86JSanszMTgDGUhl0vwRGmtkZDS3U3WcQBnvcSmjhHEcYol5SZbeHgX8Ds4EvgD9GdfWNRiL2rcNLXko4ljbRzNYCLwNDolqeB24BXkvtEz0mHZ9Zs2JasFBEpP6i4e4fAy13MOhF6kktKBGROjKzk6JzpToSTrZ+RuHU+BRQIiJ195/AMkJXYjnbH7+TRqAuPhERSSS1oEREJJF0oq6kRZcuXbywsDDuMkQkjYqLi1e4e0F9H6+AkrQoLCxk8uTJNe8oIlnDzObVvNeOKaBEJHOVlcGiRbBsGXTvDj17Qq5Wy8gWOgYlIpll8WL4859h5Eho1QoKC2H0aOjbN1wfPRpuuSWElmQ0BZSIZIaVK+G886B3b7joIsjPh0svhbvugiefhDvuCLeXlcHPfw69esGFF8JXX8VdudSTuvhEJNnc4b774JJLQtj85Cfw05/CkB2sNH/NNfDxx6EVdcst8Pe/w003wWmnVb+/JJZaUCKSXFu2wLnnhm3PPeHDD+HWW3ccTil77BFaVu+/H7oATz89tKpKS9NStjQOBZSIJNOSJXDooaH19LvfwauvhuCpi1Gj4O234Re/CK2po44KXYWSERRQIpI8ixeHcJo2DR5/HP7wB8ip55+rvDz4y1/gwQfh3XdhzBiFVIZQQIlIsnz5JRx2WLh88UX4zndqfEitnHkmPP00fPZZCKkV9V20WNJFASUiybFiBRx+eAinF16Agw5q3Oc/6qgQUjNmwJFHwtq1jfv80qgUUCKSDJs3w4knwty58NxzcOCBTfM6Rx0F//xn6D485RQNnEgwBZSIxK+iAn70I3jnnXCs6OCDm/b1xo4N5029+GIYsq5VHRJJ50GJSPz+53/C+UrXXBNaNelw7rnwxRdw7bUwdGgY6SeJohaUiMTr+efDKL0zzwwzQ6TTVVeFbsWLL4a33krva0uNFFAiEp+5c+GMM8JJuHfcAWbpff2cHHjgARgwAE49NZx7JYmhgBKReGzZAt/9bjj+9Pjj0KZNPHW0bx9ef82aMB1SWVk8dch2FFAiEo/LL4fi4tCCGTQo3lr23BPuvBPeeAOuvjreWuRrCigRSb8XXoAbb4Tzzw/HgJLgzDPhBz8Ix8PeeSfuagQw1/BKSYOioiLXiroCwNKlsNde0LUrvPcetG4dd0WV1q6FffaB8nL46KPQ/Sf1ZmbF7l5U38erBSUi6eMO55wTguBvf0tWOAHsuis8/HBYpfcnP4m7mmZPASUi6XPvvfCvf4Vzj+o6M3m67LtvmD394YfhH/+Iu5pmTV18khbq4hNmz4YRI8KS7C+9VP/ZydOhrCxMtTRrVlj8sEePuCvKSOriE5HkKy+Hs88OoXT//ckOJwhLdDz4IGzcGGac0Bf5WCT8f4mIZIWbbw4zNdx6K/TtG3c1tTNkCFx/fZi49r774q6mWVJAiUjT+uwzuOIKOOGEMJQ7k5x/Pnzzm2Gevnnz4q6m2VFAiUjTKSuDs86Ctm3jmcqooXJyQuspNfpQXX1ppYASkabzpz+Fc51uuw26d4+7mvopLIQ//xleeSWErKSNAkpEmsbHH8Pvfw8nn5y+JTSayn/8R1jo8OKLw2hESQsFlIg0vtLS0LXXvj3cfnvmde1tywzuuQdyc8PCihUVcVfULCigRKTxXXMNfPBBCKeCgriraRx9+sBNN8Gbb4bRiNLkFFAi0rimTAkr5J5+eujeyyZnnw3f/jZcdhnMnBl3NVlPASUijWfzZvjhD0OrKRtbGWZw991hDsEf/lBrRzUxBZSINJ7f/S4Mjrj3XujcOe5qmkbPnqHrctIkuO66uKvJagooEWkcb70VhpWPGwdHHx13NU3r1FPDduWVoUtTmoQmi5W00GSxWW7tWth77/Dz1KnQrl289aTDypVhJd6OHWHy5OQtHZIAmixWROJ3/vkwfz783/81j3CC0IX5wAMwfXo4P0oanQJKRBrm4YdDMP32t3DAAXFXk15HHRXm6fvf/4Vnn427mqyjLj5JC3XxZam5c8MaT3vsAW+8EZapaG42bw6LHC5eHLo3M3VKpyagLj4RiUdJSRgoAKEF1RzDCaBVq7B8/fr1cMYZYe0raRQKKBGpn8suCxPB3ncf9O8fdzXxGjYsdPO9+ir88Y9xV5M1FFAiUndPPQU33gg/+1n2zRZRX2efHU7e/cMfQlBJgymgRKRuPv88TAQ7ahTccEPc1SSHWVhWZOhQOO00WLAg7ooyngJKRGpv3bqwMm5eHjz2GLRsGXdFydK2LTzxRBg4cdJJsGlT3BVlNAWUiNRORUXowpo5Ex59NCzkJ9sbOjQMGikuhvPO0yq8DaCAEpHa+d3v4Mknw+qyhx8edzXJdvzx4VjUgw+G6Z+kXprpuFARqZO774arroJzz4ULLoi7mszwm9/AJ5/AJZdA376VQ/Kl1hRQIrJzzz8P//VfMHZsGASQ6avjpktODowfD19+GbpGe/aEgw+Ou6qMoi4+EdmxCRPge98Lk6I++ijk58ddUWZp1SoMye/fP3T7ffhh3BVlFAWUiFSvuDgsm9GzJzz3HOyyS9wVZaZOneDFF8O/35FHhsllpVYUUCKyvY8+ChOhduoEr7wCPXrEXVFm69cvnLybnw9jxmi5+FpSQInI1iZMgMMOgzZtQjj16RN3Rdlh0CB4+eUwV99BB2mhw1pQQIlIpRdfhCOOgIICePttGDAg7oqyy7BhYeXhVq3Cl4C33oq7okRTQIlIcOedcOyxMGRI+MPZr1/cFWWnIUPgnXfCsb0jjwwn9Uq1FFAizV1paVgR97zzwh/M11+Hbt3iriq79ekTvgTstx+ceSZceqmW6aiGAkqkOZszJ3Q13XZbOKH0mWegffu4q2oeunSBl14K55hdf30YlLJoUdxVJYoCSqQ5coeHHoK994aPPw4L7l13HeTmxl1Z85KfH74c3HsvTJwIe+0VJpsVQAEl0vzMng3HHAM/+EE4Afejj8LyEBKfH/84jOrr3z+sr3XyybBwYdxVxU4BJdJcrF0LV14Jw4eH4x833hiON2lW8mQYPDgM8b/66nBi9O67w7XXwsaNcVcWGwWUSLZbvz7MQD5gQJhh+7jj4NNP4cILw7pOkhwtWsDll4dJZr/5zfDzwIHw1782y6BSQIlkqzlz4OKLoXdvuOgiKCqC998Pc+r17h13dbIzAwbA00+Hc9EGD4af/SzMiP7rXzerlXoVUCLZZNmysDTGIYeEP3I33hjm05s4EV54IYSUZI4DDwzdsG+8ET7Ta64J56eNGQMPPAArV8ZdYZMy12qPkgZFRUU+efLkuMvIPps3w6RJ4Q/Y88+Hn93Dqq5nnhk2TVWUPebMCYsgPvhgGOySkwMHHBC+hBx6aPgC0rJl3FV+zcyK3b3e34oUUJIWCqgGqqiAxYthxgz47LMw8u6DD2DqVCgpCWs0jRwZji8ddxzss4/Wbcpm7qG79plnwvbRR+H2Vq1gxIjw+Y8YEb6oDBkC3bvH8v9BASUZQQFFCJmSktDq2bQpbBs2hG3tWvjqK1i9OnTbrFgBS5aEbf78cNxhy5bK5+rYMQTSyJFhEbwDDwwzj0vztGJFOF711lvhi8uUKeH/U0qrVuEYVp8+Iay6dw8nCnfuHP4vtW8flgNp1y5MEtymTXhMq1Zh4EZO/Y4GKaAkIxS1bu2TBw2Ku4ztVff/P3Vb1fvcd75VVISpalKXqa2sLEwlVFpat6ls2rQJf0S6dQt/VPr1C8PBhwwJW69eaiHJjrmH86g++yy0uufOhXnzwhedpUvDF5/Nm2v/fLm54aTi/Pww8jM3d+stJydsZmG7/nr4zncaHFAaYyrp0apVGI2URNX9oU/dVvW+1C9fdVvVX9LUL21eXuXWokX45W7RAlq3DscJ2rat/Lbavj3sumtoBXXqFPYRqS+z8MWmT58wv2J1Nm6EVavCtm5daHFt2BBu37gxBNjmzaHVX1ISvmSlvnBV/RJWUVH5xSz1ha2goFHehgJK0mPgQHj88birEJGU1JejBJ9yoGHmIiKSSDoGJWlhZuuAGXHX0Ui6ACviLqKRZMt7yZb3Adn1Xoa4+y71fbC6+CRdZjTkYGmSmNlkvZdkyZb3Adn3XhryeHXxiYhIIimgREQkkRRQki53xV1AI9J7SZ5seR+g9/I1DZIQEZFEUgtKREQSSQElIiKJpICSJmVmY81shpnNMrPL4q6nLsysj5m9ZmbTzewTM/t5dHsnM3vJzD6PLjvGXWttmVmumU0xs2ej6/3NbFL0+TxiZi3irrE2zKyDmT1mZp+Z2admtn+mfi5m9ovo/9fHZvY3M2uVKZ+Lmd1nZsvM7OMqt1X7OVhwS/SepprZyJqeXwElTcbMcoH/BY4GhgGnm9mweKuqkzLgV+4+DNgPOD+q/zLgFXffDXglup4pfg58WuX6dcCN7j4IWA2cE0tVdXcz8IK7DwVGEN5Txn0uZtYLuAAocvc9gFzgNDLnc3kAGLvNbTv6HI4Gdou2ccDtNT25Akqa0mhglrvPdvcS4O/ACTHXVGvuvtjdP4h+Xkf4I9iL8B7GR7uNB06Mp8K6MbPewDHAPdF1Aw4HHot2yYj3YmbtgUOAewHcvcTd15ChnwthwoTWZpYHtAEWkyGfi7u/Caza5uYdfQ4nAA96MBHoYGY9dvb8DQqoTO6+kbToBSyocn1hdFvGMbNCYB9gEtDN3RdHdy0BusVUVl3dBFwCVETXOwNr3L0sup4pn09/YDlwf9RdeY+ZtSUDPxd3XwT8CZhPCKavgGIy83NJ2dHnUOe/B/UOqCzovhGpFTNrBzwOXOjua6ve5+E8jcSfq2FmxwLL3L047loaQR4wErjd3fcBNrBNd14GfS4dCS2L/kBPoC3bd5llrIZ+Dg1pQWV0942kxSKgT5XrvaPbMoaZ5RPC6SF3fyK6eWmqayK6XBZXfXVwIHC8mc0l/K4eTjiO0yHqWoLM+XwWAgvdfVJ0/TFCYGXi53IEMMfdl7t7KfAE4bPKxM8lZUefQ53/HtT7RF0z+y4w1t3Pja6fCezr7j/d0WO6dOnihYWF9Xo9Sb65c8O6Z3vuuf19xcXFKwn/Gb/v7p+kubR6iY7RjAdWufuFVW6/4Qj77kXxVSbZ7GV/7ALCMbbH3f3vZnYHMNXdb4u5tGpF3d/PRoM8MLMbgJXufm106KeTu19iZscAPwW+DewL3OLuo3f23E0+m7mZjSOM2KBv375MntygyW0lwU49FaZOheo+4qib7NFMCafIgcCZwDQz+zC67QrgWkABJU3lLuBfwN/N7I/AFKIBIUljZn8DDgO6mNlC4PeE349HzewcYB5wSrT7c4RwmgVsBH5U0/M3JKBq1Vxz97uI5mMqKipKfJ+w1N/69WGBzh342N2vSmM5DebubwPVrAcPR+Z8L83V1MAsLLUtGc/dtwCzCYdREs3dT9/BXWOq2deB8+vy/A05BvU+sFt0QlkLwtj9pxvwfJLhVq2CjhlxaqSIZIJ6B1Q0BPKnwIuE80MyrftGGtnChdC7d9xVNFPuoRVlBjm5lT+LZLAGHYNy9+cI/YrSzK1aFQJqyJC4K2nGUl18Xh5CCrC8XLy8PNxeUR5TYSL1o5kkpFG89lq43H//eOsQkezR5KP4pHm45x7o2RMOOijuSgT4urXkFeVYfphn1Nq2wTdvCbeXlsRWmkhtKaCkwZ56Cl54Aa6+GvL0PypxUmHkpSXkRMMsczp1wDduAqBi3brYahPZGXXxSYNMnAg/+hGMGAG/+lXc1YhINlFASb24w8MPw5gx0KkT/POf0CKRK9ZIVRUbN1KxcSPlS5dheXlYXh65gweSW1BAbkFB3OWJbEUdMlJnU6fCBRfAG2/A6NGhi69797irkroqX706/LB6NXl9ovMDRu9J7qr1AFTMW6RjVRIrtaCkVkpK4B//gCOPDN1506bBHXfAhAkKJxFpGmpByQ5t2QJvvQXPPQcPPQTLlkHfvvDf/w0/+Ql07hx3hdJYyhYsBMCWLYdBhQCUHrIneRtKAcj9YjHly5fHVZ40Uwoo+Zo7zJoFL74YRuW99hps3BiOLX372zBuHBx1FOTmxl2pNBXfsoXyT2YA0GplN7YMDevJbTx4AFbRH4B2M1ZT/unnsdUozYcCqhlbuhTefz9s770XLleuDPftthuccw6MHQuHHgpt28Zbq4g0PwqoZqCkJLSMpk+HTz+Fjz4KYTR/frg/JweGD4cTTgiDHo44AgYOjLdmiV/ZkqXkLlkKwC7Dh7B6RJgJeNmBXeCALgC0nxMGUbSY+CkVGzfGU6hkLQVUFtm0CWbOrAyi6dPD9vnnUFZWud+AAWFKogsuCIG0zz7Qrl18dUvylX8yg12nh8lny745khV7tgRgyb7hkn33pu2XYS7ALu8uo3zmF7HUKdlFAZVB3EO33OzZYfvii61/Xry4ct/c3NAKGjYMTjoJdt89/DxkiLrrRCQzKKASxD2MlJs/f+ttzpzKMKrai2IWlrcYMACOPjpc7rZbCKLddoOWLeN7L5KFotnS814tpvfMMHhi5aFhzdLVuxsriioAWLVXAa2WdwOg2/tbyHulOIZiJRsooNJo40ZYsGD7AEptCxaEod1VtWkTgmfAgMpjQwMGhMt+/aBVq3jeizRvZQvD4tntH4ou99uLRYeGfuIN/cuo6BdO9p27Rz4cux8AnT80ujwbRgiWr1yV7pIlAymgGsnatWE9pG23RYsqf161ze+kWZgBvG9fGDUqdMX17bv11rGj1p0TkeZJAVUDd1i9uvrwqbpVNyF0166hC66wMCxD0avX1uHTs6fmr5MsMXEqvSaGHzd+Z1++PCS0ptoWfsXgfgsAWDG8HTO/1ReAltOG0ve5MNVSxUefpr9eyQjNOqAqKmD58h2HTqr1s2nT1o/LyQnT+/TpE473HHVUCKLevUMI9e4dwkfHgKQ5avPEJIa8GaYZWXzKED46KPwijOo3n6O6hzD6avfWvHPogLDPZ/vR77kwzDT/35NjqFiSKqsDqqQkHNeZNy8c45k3b+ufqzvmk5dXGTIjR8Lxx1eGT2rr3l3rHomINLWM/jObGvU2c2bllgqhefNgyZKvBx59rUeP0L02cmQ45tOnT9hS4dO1a2ghiUj9la8IU5J0vW0C3d4dDsCHJw1h+QGh6++Y7h9z+aDnACgdmMdjo4sAmPD9UfT4V/iz1O4fk9JdtiRMRgTUunXhZNOZM2HGjK0Dae3ayv3y88PItr59wxQ9qZ9Tl336qNtNJN18yicA9JsCJd8KQXTrCUdw8D6fAfCf3V7n6t7PAlDaC+4acRAAz542nLbP7ApAxwfeTXfZkgCJC6jUhKXvvBO2CRPCbAgpZiFsBg+GM88Ml6mtXz9NZCoiki3Mt+0D23YHsz7Ag0A3wIG73P1mM+sEPAIUAnOBU9x99c6eq6ioyCdP3v4gaEUFjB8fFr6bMCEMXADo0CFMybP//mGuuMGDw/k/rVvX9W1K3Mys2N2L4q6jsRyZ872d/+JIjdadFs6PWnb8Zs4Y/j4AP+v0HrvkhKGtpV7OXWuGAfDQnCJy/hkGXnS6L7tbUy9V/EMnlkRqE1A9gB7u/oGZ7QIUAycCZwOr3P1aM7sM6Ojul+7suaoLqFmz4Nxzw+qsAwbAwQfDgQfCAQeE6Xl0PCg7KKBkZ1aftT8AW05aw1mDwrGnX3aa/fX9W7yUW1YPBWD8zP1o83TU9Td+4vYHmjOcAqpSjV187r4YWBz9vM7MPgV6AScAh0W7jQdeB3YaUNvasAGKiuCrr+CWW+CnP9VJqSIiEtTYgtpqZ7NC4E1gD2C+u3eIbjdgder6No8ZB4wD6Nu376h58+Z9fZ87XHwx/OUv4WTWyy8PLaghQxRU2UYtKKmtlf8RWlOlx6zh+wNDj8v5HafR2kLX3+qKTdy6ajQAD33yDbo9Geb7avfoxBiqbXxqQVWqdUCZWTvgDeAqd3/CzNZUDSQzW+3uHXf2HDs6BvX226Gbb0aYpotOnUIXX6qr7xvf0HGnTKeAkvpYfXYIqw3HreWEAdMAOK/zBMqjf/0F5e24Z+khALw1ZSiFT4UJa1u8mLkn/CqgKtVqFJ+Z5QOPAw+5+xPRzUvNrIe7L46OUy2rbxEHHRTWL5oxIwySSI3gezaMPCU3NxyfqjpiL7X16qXWlohINqrNIAkjHGNa5e4XVrn9BmBllUESndz9kp09145aUDuyYgW8+y5MmlR5DtTnn2899VCbNmFpiVRg9e+/9blPmu07GdSCkoZaf0oY9bf4mBIOGjwLgO90+YBWVgrA9M29ePrLvQBY8m5P+v9zDQAVH06v5tmSSy2oSrUJqIOAt4BpQEV08xXAJOBRoC8wjzDMfKdz6Nc1oKpTURHmyKt6sm5qmzMHysu33r9r1xBW2560m7pNs4WnhwJKGlPZmFEAzD0mn17Dw7L0+xbMpWNeWDBt6rpeFM8LE9O2e7sNAD0e+YyKtWEZEC8tSXfJtaaAqlSbUXxvAzv6BxvTuOXULCencnqiMdu8emlpCK9t59ybNw+mTQtdhps3b/2Ytm23n2tv261zZ4WYiEi6JW4miYbIzw+jAQsLq7/fPZwEnAqu1GVq1vJXXoEvvwyttKpatqw5xDSHn0j6pFbpHfQK5OwVzo967tj9Kd1zAwB9ClYzvNdiAGYf1QmAGbsNptcb4Ze79ZPvpbtkqYesCqiamIUg6do1nH9VnbIyWLp0+8UGU9uECeGytHTrx6VmQU/NhF7d1qOHZkEXEakt/bncRtWg2ZGKijCAY0frSE2ZAs88U/06UqnZ1He06ZiYSN1UTA2TzvaeCrmdQ2tpxbFDmDMqtJasUzje1LLPehZ8KxyPajf4AHq+Hq0y+t60NFcstVWnE3UbqjEGSWSK1Eq8VVthCxaEbf78yq1km2O1bdtuH1r9+4dh9gMHQkFBZgaYBklIXEqPCAMqln6jJVs6hY+tonUF+WtCn3yn6U6nN8Oqv2ULF8VTZBUaJFFJLagmYhZOOO7UCfbcs/p9Uiv6Vg2sqtuUKWG9q6ratQthldoGDqy87NdPS8iLSPZQQMUoJwe6dQvbN75R/T6bNoWBHF98AbNnV15+/jm8+OLW3Yg5OSGkhg0L2+67V17uumt63pNI0uS/HAZU9H4ZcgcPBGDF/l3Z0DM0VNYWGpu69AOg88fdyXu1OJ5CZTsKqIRr3RqGDg3bttzDqsGp4Prii3A+2PTp8PLLWy9n36tXZXANHx4CcfjwMPJRpLkon/kFAB1nfkHntm0BKNlvKF8Vhq6Hr/q3IP/74YTgjh+uonz6zHgKFUABldHMwqCLHj3CvIVVlZeHE5enT996u/tu2BjUbs6jAAAIEklEQVTOZaRVKxg5MoTV6NHhctCgzDzGJSLZR4MkmpmKitDiev/9sL33HnzwQWVXYceO4QTosWPDtrPRjHWhQRKSCXKHDWb94DAHtucYLVeF80lafLqQ8qX1nm60TjRIopJaUM1MTk5oJQ0aBKefHm4rK4NPPgmBNWEC/Pvf8Nhj4b499wxBdfLJoZWl1pVks/LpM2kdTd2X260r5f27A7Bljz607BIWa6iYOSfRUyVlE819IOTlwYgRYcmT++4LQ+GnTYMbbggnNd90E+y3H+y9N9x2W1hgUkSkqSmgZDtmsMcecNFFYbDFihVw550hyM4/PxzzuuIKWL8+7kpFmk750mUwcSpMnEr+m9OgpBRKSvFRQ8nr1ZO8Xj3jLjHrKaCkRrvuCuPGQXExTJ4MJ50E11wTRhamugJFspmXllD++WzKP58NE6fiZWV4WRm5gweS26E9uR3ax11iVlJASZ2MGgUPPRQWlOzaFb73PbjqqjDkXUSkMWmQhNTLAQeEhSR//GP4zW9CQP3mN3FXJZIeX4/oW7qM3Ogs+NxuXfGv1gJQse26PlIvCiipt/x8GD8+nHN15ZVhtN+OZokXyVbla0MosXYtOdES3jm77IJH5254WVlcpWU8dfFJg+TkhJF9BQXw61/HXY2IZBMFlDRYhw5w3nnw0kthkluR5qpi8+awrVv39W2W3yIMjdVJhHWmgJJGcfLJ4TjU66/HXYlIMqRG+nlpCVhO2HJy4y4royigpFEMGxbWsirWRNAi0kg0SEIaRU5OWFhx3ry4KxFJoIryyp9TXX06N6NGtW5BmVmumU0xs2ej6/3NbJKZzTKzR8xMS+U1c127hgUYRWQn3BVOtVSXLr6fA59WuX4dcKO7DwJWA+c0ZmGSedq3h9SIWxGRhqpVQJlZb+AY4J7ougGHA6mJbsYDJzZFgZI5WrfeeoVfEZGGqG0L6ibgEqAiut4ZWOPuqTPQFgLVrhxkZuPMbLKZTV6u/p+s1qIFlJbGXYWIZIsaA8rMjgWWuXu9xme5+13uXuTuRQUFBfV5CskQublhVgkRkcZQm1F8BwLHm9m3gVbArsDNQAczy4taUb2BRU1XpmSCnBwFlIg0nhpbUO5+ubv3dvdC4DTgVXc/A3gN+G6021nAU01WpWQEnSgvIo2pISfqXgr80sxmEY5J3ds4JUkm0+hZEWksdTpR191fB16Pfp4NjG78kiRTmSmgRKTxaKojaTTq4hORxqSAEhGRRFJAiYhIIimgREQkkRRQIiKSSAooERFJJAWUiIgkkgJKREQSSQElIiKJpIASEZFEUkCJiEgiKaBERCSRFFAiIpJICigREUkkBZSIiCSSAkpERBJJASUiIomkgBIRkURSQImISCIpoEREJJEUUCIikki1Cigz62Bmj5nZZ2b2qZntb2adzOwlM/s8uuzY1MWKiEjzUdsW1M3AC+4+FBgBfApcBrzi7rsBr0TXRUREGkWNAWVm7YFDgHsB3L3E3dcAJwDjo93GAyc2VZEiItL81KYF1R9YDtxvZlPM7B4zawt0c/fF0T5LgG7VPdjMxpnZZDObvHz58sapWkREsl5tAioPGAnc7u77ABvYpjvP3R3w6h7s7ne5e5G7FxUUFDS0XhERaSZqE1ALgYXuPim6/hghsJaaWQ+A6HJZ05QoIiLNUY0B5e5LgAVmNiS6aQwwHXgaOCu67SzgqSapUEREmqW8Wu73M+AhM2sBzAZ+RAi3R83sHGAecErTlCgiIs1RrQLK3T8Eiqq5a0zjliMiIhJoJgkREUkkBZSIiCSSAkpERBJJASUiIomkgBIRkURSQImISCIpoEREJJEUUCIikkgKKBERSSQFlIiIJJICSkREEkkBJSIiiaSAEhGRRFJAiYhIIimgREQkkRRQIiKSSAooERFJJAWUiIgkkgJKREQSSQElIiKJpIASEZFEqlVAmdkvzOwTM/vYzP5mZq3MrL+ZTTKzWWb2iJm1aOpiRUSk+agxoMysF3ABUOTuewC5wGnAdcCN7j4IWA2c05SFiohI81LbLr48oLWZ5QFtgMXA4cBj0f3jgRMbvzwREWmuagwod18E/AmYTwimr4BiYI27l0W7LQR6Vfd4MxtnZpPNbPLy5csbp2oREcl6teni6wicAPQHegJtgbG1fQF3v8vdi9y9qKCgoN6FiohI81KbLr4jgDnuvtzdS4EngAOBDlGXH0BvYFET1SgiIs1QbQJqPrCfmbUxMwPGANOB14DvRvucBTzVNCWKiEhzVJtjUJMIgyE+AKZFj7kLuBT4pZnNAjoD9zZhnSIi0szk1bwLuPvvgd9vc/NsYHSjVyQiIoJmkhARkYRSQImISCIpoEREJJEUUCIikkgKKBERSSQFlIiIJJICSkREEkkBJSIiiaSAEhGRRFJAiYhIIimgREQkkRRQIiKSSAooERFJJAWUiIgkkgJKREQSSQElIiKJpIASEZFEUkCJiEgiKaBERCSRFFAiIpJICigREUkkBZSIiCSSAkpERBJJASWNpl8/2HvvuKsQkWxh7p6+FzNbDsxL2wtKkvRz94K4ixCRzJHWgBIREaktdfGJiEgiKaBERCSRFFAiIpJICigREUkkBZSIiCSSAkpERBJJASUiIomkgBIRkURSQImISCIpoEREJJEUUCIikkgKKBERSSQFlIiIJJICSkREEkkBJSIiiaSAEhGRRFJAiYhIIimgREQkkRRQIiKSSAooERFJJAWUiIgkkgJKREQS6f8Dgok+PewAtxcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% Example with Frobenius norm regularization\n",
    "\n",
    "\n",
    "def f(G):\n",
    "    return 0.5 * np.sum(G**2)\n",
    "\n",
    "\n",
    "def df(G):\n",
    "    return G\n",
    "\n",
    "\n",
    "reg = 1e-1\n",
    "\n",
    "Gl2 = ot.optim.cg(a, b, M, reg, f, df, verbose=True)\n",
    "\n",
    "pl.figure(3)\n",
    "ot.plot.plot1D_mat(a, b, Gl2, 'OT matrix Frob. reg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD with entropic regularization\n",
    "--------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "    0|1.692289e-01|0.000000e+00\n",
      "    1|1.617643e-01|-4.614437e-02\n",
      "    2|1.612639e-01|-3.102965e-03\n",
      "    3|1.611291e-01|-8.371098e-04\n",
      "    4|1.610468e-01|-5.110558e-04\n",
      "    5|1.610198e-01|-1.672927e-04\n",
      "    6|1.610130e-01|-4.232417e-05\n",
      "    7|1.610090e-01|-2.513455e-05\n",
      "    8|1.610002e-01|-5.443507e-05\n",
      "    9|1.609996e-01|-3.657071e-06\n",
      "   10|1.609948e-01|-2.998735e-05\n",
      "   11|1.609695e-01|-1.569217e-04\n",
      "   12|1.609533e-01|-1.010779e-04\n",
      "   13|1.609520e-01|-8.043897e-06\n",
      "   14|1.609465e-01|-3.415246e-05\n",
      "   15|1.609386e-01|-4.898605e-05\n",
      "   16|1.609324e-01|-3.837052e-05\n",
      "   17|1.609298e-01|-1.617826e-05\n",
      "   18|1.609184e-01|-7.080015e-05\n",
      "   19|1.609083e-01|-6.273206e-05\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   20|1.608988e-01|-5.940805e-05\n",
      "   21|1.608853e-01|-8.380030e-05\n",
      "   22|1.608844e-01|-5.185045e-06\n",
      "   23|1.608824e-01|-1.279113e-05\n",
      "   24|1.608819e-01|-3.156821e-06\n",
      "   25|1.608783e-01|-2.205746e-05\n",
      "   26|1.608764e-01|-1.189894e-05\n",
      "   27|1.608755e-01|-5.474607e-06\n",
      "   28|1.608737e-01|-1.144227e-05\n",
      "   29|1.608676e-01|-3.775335e-05\n",
      "   30|1.608638e-01|-2.348020e-05\n",
      "   31|1.608627e-01|-6.863136e-06\n",
      "   32|1.608529e-01|-6.110230e-05\n",
      "   33|1.608487e-01|-2.641106e-05\n",
      "   34|1.608409e-01|-4.823638e-05\n",
      "   35|1.608373e-01|-2.256641e-05\n",
      "   36|1.608338e-01|-2.132444e-05\n",
      "   37|1.608310e-01|-1.786649e-05\n",
      "   38|1.608260e-01|-3.103848e-05\n",
      "   39|1.608206e-01|-3.321265e-05\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   40|1.608201e-01|-3.054747e-06\n",
      "   41|1.608195e-01|-4.198335e-06\n",
      "   42|1.608193e-01|-8.458736e-07\n",
      "   43|1.608159e-01|-2.153759e-05\n",
      "   44|1.608115e-01|-2.738314e-05\n",
      "   45|1.608108e-01|-3.960032e-06\n",
      "   46|1.608081e-01|-1.675447e-05\n",
      "   47|1.608072e-01|-5.976340e-06\n",
      "   48|1.608046e-01|-1.604130e-05\n",
      "   49|1.608020e-01|-1.617036e-05\n",
      "   50|1.608014e-01|-3.957795e-06\n",
      "   51|1.608011e-01|-1.292411e-06\n",
      "   52|1.607998e-01|-8.431795e-06\n",
      "   53|1.607964e-01|-2.127054e-05\n",
      "   54|1.607947e-01|-1.021878e-05\n",
      "   55|1.607947e-01|-3.560621e-07\n",
      "   56|1.607900e-01|-2.929781e-05\n",
      "   57|1.607890e-01|-5.740229e-06\n",
      "   58|1.607858e-01|-2.039550e-05\n",
      "   59|1.607836e-01|-1.319545e-05\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   60|1.607826e-01|-6.378947e-06\n",
      "   61|1.607808e-01|-1.145102e-05\n",
      "   62|1.607776e-01|-1.941743e-05\n",
      "   63|1.607743e-01|-2.087422e-05\n",
      "   64|1.607741e-01|-1.310249e-06\n",
      "   65|1.607738e-01|-1.682752e-06\n",
      "   66|1.607691e-01|-2.913936e-05\n",
      "   67|1.607671e-01|-1.288855e-05\n",
      "   68|1.607654e-01|-1.002448e-05\n",
      "   69|1.607641e-01|-8.209492e-06\n",
      "   70|1.607632e-01|-5.588467e-06\n",
      "   71|1.607619e-01|-8.050388e-06\n",
      "   72|1.607618e-01|-9.417493e-07\n",
      "   73|1.607598e-01|-1.210509e-05\n",
      "   74|1.607591e-01|-4.392914e-06\n",
      "   75|1.607579e-01|-7.759587e-06\n",
      "   76|1.607574e-01|-2.760280e-06\n",
      "   77|1.607556e-01|-1.146469e-05\n",
      "   78|1.607550e-01|-3.689456e-06\n",
      "   79|1.607550e-01|-4.065631e-08\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   80|1.607539e-01|-6.555681e-06\n",
      "   81|1.607528e-01|-7.177470e-06\n",
      "   82|1.607527e-01|-5.306068e-07\n",
      "   83|1.607514e-01|-7.816045e-06\n",
      "   84|1.607511e-01|-2.301970e-06\n",
      "   85|1.607504e-01|-4.281072e-06\n",
      "   86|1.607503e-01|-7.821886e-07\n",
      "   87|1.607480e-01|-1.403013e-05\n",
      "   88|1.607480e-01|-1.169298e-08\n",
      "   89|1.607473e-01|-4.235982e-06\n",
      "   90|1.607470e-01|-1.717105e-06\n",
      "   91|1.607470e-01|-6.148402e-09\n",
      "   92|1.607462e-01|-5.396481e-06\n",
      "   93|1.607461e-01|-5.194954e-07\n",
      "   94|1.607450e-01|-6.525707e-06\n",
      "   95|1.607442e-01|-5.332060e-06\n",
      "   96|1.607439e-01|-1.682093e-06\n",
      "   97|1.607437e-01|-1.594796e-06\n",
      "   98|1.607435e-01|-7.923812e-07\n",
      "   99|1.607420e-01|-9.738552e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  100|1.607419e-01|-1.022448e-07\n",
      "  101|1.607419e-01|-4.865999e-07\n",
      "  102|1.607418e-01|-7.092012e-07\n",
      "  103|1.607408e-01|-5.861815e-06\n",
      "  104|1.607402e-01|-3.953266e-06\n",
      "  105|1.607395e-01|-3.969572e-06\n",
      "  106|1.607390e-01|-3.612075e-06\n",
      "  107|1.607377e-01|-7.683735e-06\n",
      "  108|1.607365e-01|-7.777599e-06\n",
      "  109|1.607364e-01|-2.335096e-07\n",
      "  110|1.607364e-01|-4.562036e-07\n",
      "  111|1.607360e-01|-2.089538e-06\n",
      "  112|1.607356e-01|-2.755355e-06\n",
      "  113|1.607349e-01|-4.501960e-06\n",
      "  114|1.607347e-01|-1.160544e-06\n",
      "  115|1.607346e-01|-6.289450e-07\n",
      "  116|1.607345e-01|-2.092146e-07\n",
      "  117|1.607336e-01|-5.990866e-06\n",
      "  118|1.607330e-01|-3.348498e-06\n",
      "  119|1.607328e-01|-1.256222e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  120|1.607320e-01|-5.418353e-06\n",
      "  121|1.607318e-01|-8.296189e-07\n",
      "  122|1.607311e-01|-4.381608e-06\n",
      "  123|1.607310e-01|-8.913901e-07\n",
      "  124|1.607309e-01|-3.808821e-07\n",
      "  125|1.607302e-01|-4.608994e-06\n",
      "  126|1.607294e-01|-5.063777e-06\n",
      "  127|1.607290e-01|-2.532835e-06\n",
      "  128|1.607285e-01|-2.870049e-06\n",
      "  129|1.607284e-01|-4.892812e-07\n",
      "  130|1.607281e-01|-1.760452e-06\n",
      "  131|1.607279e-01|-1.727139e-06\n",
      "  132|1.607275e-01|-2.220706e-06\n",
      "  133|1.607271e-01|-2.516930e-06\n",
      "  134|1.607269e-01|-1.201434e-06\n",
      "  135|1.607269e-01|-2.183459e-09\n",
      "  136|1.607262e-01|-4.223011e-06\n",
      "  137|1.607258e-01|-2.530202e-06\n",
      "  138|1.607258e-01|-1.857260e-07\n",
      "  139|1.607256e-01|-1.401957e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  140|1.607250e-01|-3.242751e-06\n",
      "  141|1.607247e-01|-2.308071e-06\n",
      "  142|1.607247e-01|-4.730700e-08\n",
      "  143|1.607246e-01|-4.240229e-07\n",
      "  144|1.607242e-01|-2.484810e-06\n",
      "  145|1.607238e-01|-2.539206e-06\n",
      "  146|1.607234e-01|-2.535574e-06\n",
      "  147|1.607231e-01|-1.954802e-06\n",
      "  148|1.607228e-01|-1.765447e-06\n",
      "  149|1.607228e-01|-1.620007e-08\n",
      "  150|1.607222e-01|-3.615783e-06\n",
      "  151|1.607222e-01|-8.668516e-08\n",
      "  152|1.607215e-01|-4.000673e-06\n",
      "  153|1.607213e-01|-1.774103e-06\n",
      "  154|1.607213e-01|-6.328834e-09\n",
      "  155|1.607209e-01|-2.418783e-06\n",
      "  156|1.607208e-01|-2.848492e-07\n",
      "  157|1.607207e-01|-8.836043e-07\n",
      "  158|1.607205e-01|-1.192836e-06\n",
      "  159|1.607202e-01|-1.638022e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  160|1.607202e-01|-3.670914e-08\n",
      "  161|1.607197e-01|-3.153709e-06\n",
      "  162|1.607197e-01|-2.419565e-09\n",
      "  163|1.607194e-01|-2.136882e-06\n",
      "  164|1.607194e-01|-1.173754e-09\n",
      "  165|1.607192e-01|-8.169238e-07\n",
      "  166|1.607191e-01|-9.218755e-07\n",
      "  167|1.607189e-01|-9.459255e-07\n",
      "  168|1.607187e-01|-1.294835e-06\n",
      "  169|1.607186e-01|-5.797668e-07\n",
      "  170|1.607186e-01|-4.706272e-08\n",
      "  171|1.607183e-01|-1.753383e-06\n",
      "  172|1.607183e-01|-1.681573e-07\n",
      "  173|1.607183e-01|-2.563971e-10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XeYldW5/vHvM41eBek4IpYgqEES7EeDJtYfnmiMRo94TMTeYmJNIclRozGWGDUxEmOPirHGEjUajQoKYhRBAVEUlCIOHZm2fn887zszDAMzTFt777k/17Wvze4PG+aeZ6+93rUshICIiLS+vNgFiIi0VQpgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASzSwsxsXzN7P3YdknkUwJLxzOwkM3vHzNaa2SIzu8XMuie3/cHMVienUjMrq3H5qVaoLZjZ0M3dJ4Twcghhx0Y+/0dmtq7G32m1mf2+gY990cx+0JjXldahAJaMZmYXAFcBPwa6AXsA2wDPmllRCOG0EELnEEJn4Arg/vRyCOGQeJU7Mytohqc5osbfqXMI4axmeM7mqi1jXy8bKIAlY5lZV+AXwNkhhKdDCGUhhI+AY4Bi4IRGPOf+ZrbAzC40syVm9pmZHWlmh5rZbDP7wswurXH/r5vZa2a2PLnv782sKLntpeRu/0k60+/WeP6LzGwRcHt6XfKY7ZLXGJlc7m9mS81s/0b8XU4ys3+b2TVmVmJmH5rZIcltlwP7Ar+v2TUnHfuZZjYHmJNct5eZvWFmK5LzvWq8xotmdqWZvW5mK83sUTPr2cD6JpjZJDO728xWAieZWZ6ZXWxmH5jZMjN7oObzmdmJZjY/ue2nySeAA7f0vckWCmDJZHsB7YG/1bwyhLAaeBI4qJHP2zd53gHAz4A/4WG+Ox5aPzWzbZP7VgDnA72APYExwBlJHfsl99k16Uzvr/H8PfFOfXyt2j8ALgLuNrOOwO3AHSGEFxv5dxkNvJ/UdzUw0cwshHAZ8DJwVh1d85HJ44Yl4fd34HfAVsC1wN/NbKsa9z8ROBnoB5Qn922oscAkoDtwD3B28vr/BfQHSoCbAMxsGHAzcHzyWt3wf6OcpQCWTNYL+DyEUF7HbZ8ltzdGGXB5CKEM+GvyPDeEEFaFEN4FZgK7AoQQpoUQJocQypPu+494eGxOJfDzEML6EMK62jeGEP4EzAWm4EFzWT3P90jSgaenU2rcNj+E8KcQQgVwR/J8fep5vitDCF8ktR0GzAkh3JX8He8D3gOOqHH/u0IIM0IIa4CfAseYWX49r5F6LYTwSAihMnm904DLQggLQgjrgQnA0cnwxNHA4yGEf4cQSvFfjjm9WI3GZCSTfQ70MrOCOkK4X3J7YyxLAgsgDcjFNW5fB3QGMLMd8K5wFNAR/5mZVs/zLw0hfFnPff4EPAaMT4Joc44MITy3idsWpX8IIaw1M9LaN+OTGn/uD8yvdft8Nuw8P6l1WyH+S6vme9aQ1wL/VPCwmVXWuK4C/6XRv+b9k7/Psga8RtZSByyZ7DVgPfDtmleaWWfgEOD5VqjhFrwj3D6E0BW4FLB6HrPZri2p/3pgIjChoWOqjbCpOmpe/ykeijUNBhbWuDyo1m1lNPyXX+0aPgEOCSF0r3FqH0JYiH+qGZje0cw64MMiOUsBLBkrhLAC/xLuRjM72MwKzawYeABYANzVCmV0AVYCq81sJ+D0WrcvBoZs4XPeAEwNIfwAH3/9Q5OrrFtDansS2MHMvmdmBWb2XWAY8ESN+5xgZsOSMetfApNqfILYUn8ALjezbQDMrLeZjU1umwQckXwpWIQPT9T3yy6rKYAlo4UQrsa7zmvwIJyCd1FjGvDRvTn8CPgesAofNri/1u0TgDuSsdlj6nuyJGwOpjrIfwiMNLPjN/Owx2vNA364gbXfgI+vlphZnV+chRCWAYcDFwDLgAuBw0MINTvcu4C/4MMd7YFzavx9VpvZvg2sJ63pMeAfZrYKmIx/IUgy/n42Pi7/GbAaWIJ/CspJpgXZRWRTzOxF4O4Qwm0RXrszsBwf/vmwtV+/NagDFpGMYWZHmFlHM+uEf+p5B/goblUtRwEsIplkLP7F4KfA9sCxIYc/pmsIQkQkEnXAIiKR6EAMAaBXr16huLg4dhkiWWXatGmfhxB6N/bxCmABoLi4mKlTp8YuQySrmFntowi3iIYgREQiUQCLtCUrVsCrr8L8+aAv4KNTAIvkuhDgwQdh+HDo3h323huKi2HrreHii2HVqtgVtlkKYJFctmgRjBkDxxwD+flw+eXw6KNw881wwAFw1VWw447wzDOxK22T9CWcSK5asMDDd8ECuOkmOPVUD+HU6afDlClwyilwxBFw331w1FHx6m2D1AGL5KKFC2Hffb0D/sc/4IwzNgzf1OjR8PLL8LWveZf84IOtX2sbpgAWyTVlZR6mn38Ozz/vY76b062bD0HsuSeMGwczZrROnaIAFsk5F17oMx0mToRRoxr2mM6dYdIkD+OjjoKVK1u2RgEUwCK55Ykn4Prr4ZxzvAveEn37wv33wwcfwJlntkx9sgEFsEiuWL3ax3qHD4ff/KZxz7HffnDppXD33fDss81bn2xEASySKyZMgE8+gT/+EYqKGv88l14K22/vYf5lfXuLSlMogEVywdtv+9DDKafAXns17bnat/d5wnPnwpVXNk99UicFsEguuPBC6NoVfv3r5nm+Aw+E737XhzI+/bR5nlM2ogAWyXYvvODTyC69FHo24w73V1wB5eXwy18233PKBhTAItksBF/PYeBAOOus5n3uIUP86LnbboPZs5v3uQVQAItkt4cfhtdfh1/8wsdum9tPfuLP+9OfNv9ziwJYJGuFAP/3fz5j4cQTW+Y1+vSB887zQ5RnzWqZ12jDFMAi2eqpp2D6dLjkEihowXW1zjsPOnRovi/4pIoCWCQbpd3v4MFwwgkt+1q9esH48XDPPTBvXsu+VhujABbJRv/6F7z2mk8/Kyxs+df70Y98NbWrr27512pDFMAi2ejqq31Hi5NPbp3XGzDAV0r7y19gyZLWec02QAEskm1mzvTx37PO8rHZ1vLDH8L69X6UnDQLBbBItrn+ep8adtpprfu6O+0Ehx3mAbxuXeu+do5SAItkkyVL4M47fdpZ796t//oXXABLl/oXctJkCmCRbPKHP/gwwPnnx3n9/feH3XaDa6/VtvbNQAEski1KS+GWW+CQQ3w4IAYzD/9Zs3y7I2kSBbBItpg0yTfZPPvsuHV897s+/HHjjXHryAEKYJFsceONftjxt74Vt4527fzAjMcfhw8/jFtLllMAi2SDqVNh8mSfepaXAT+2p5/udWhKWpNkwL+kiNTrxht95+KTTopdiRswwHdPvu02WLs2djVZSwEskuk+/9x3Kz7xRN/1IlOcdRYsXw733Re7kqylABbJdBMn+tSzM86IXcmG9tkHRoyAm27SlLRGUgCLZLKKCp96tv/+sPPOsavZkJn/Upg+3cenZYspgEUy2ZNPwvz5cOaZsSup2wkn+LDI738fu5KspAAWyWQ33QT9+8PYsbErqVvnzr5K2oMPwuLFsavJOgpgkUw1Z47vdnzqqa2z5m9jnXEGlJX5jAjZIgpgkUx1yy2+1dApp8SuZPN22gnGjIE//tG3sZcGUwCLZKK1a+H2232ubb9+saup35lnwiefwBNPxK4kqyiARTLRvff6HNtM/fKttiOOgEGDfMxaGkwBLJJpQvAgGzHC59pmg4ICXyD+uefgvfdiV5M1FMAimeaVV+Ctt/xIM7PY1TTcD34ARUWakrYFFMAimebGG6F7dzj++NiVbJmtt4Zjj4U77oCVK2NXkxUUwCKZZOFCeOgh+P73oVOn2NVsubPPhtWrffdkqZcCWCST/OEPUFmZees+NNSoUTB6tA9DVFbGribjKYBFMsWXX/pc2sMPhyFDYlfTeGef7QeRPP107EoyngJYJFPcc4/vOHzeebEraZrvfMfnLl93XexKMp4CWCQThOCBtcsucMABsatpmqIin8Hx3HPwzjuxq8loCmCRTPDcc/Duu77jcDZNPduUU0+FDh3ghhtiV5LRFMAimeC666BPHzjuuNiVNI+ttvIdPO6+G5YsiV1NxlIAi8T2zjvw1FN+2HG7drGraT7nnw+lpdq+fjMUwCKxXX21z/nNlnUfGmrHHX0d45tu8rnBshEFsEhM8+f7ppannAI9e8aupvlddBGUlMCf/hS7koykABaJ6dpr/Uu3H/4wdiUtY4894L/+y/+epaWxq8k4CmCRWJYs8c7w+ON9KcdcddFFsGAB3HVX7EoyjgJYJJbf/Ma3m7/kktiVtKyDD4bdd4fLL/eti6SKAlgkhiVL4Oab4Xvf8y+rcpkZTJgAH36oLrgWBbBIDNdc42s//OQnsStpHYcdpi64Dgpgkda2aJFPzTruuNzvflNpFzxvnpaqrEEBLNLafvELnxEwYULsSlrXYYf5rIgJE3zTUVEAi7Sq2bN95sOpp8LQobGraV1mftDJp59qjYiEAlikNV16qS9S87Ofxa4kjn339R2Uf/1rWLYsdjXRKYBFWstLL/l2Qz/6ke+f1lZdeaUfmvzzn8euJDoFsEhrKC/3NXK32QZ+/OPY1cS1886+5dItt/juz22YAlikNdx0k696dt110LFj7Gri+9WvfMnKM89s03vHKYBFWtqnn/qY7ze/CUceGbuazNC9O1x1Fbz6qm9j30YpgEVaUgg+46G01HcKzoXdLprLuHH+pdz558PChbGriUIBLNKS7r4bnngCrrgCtt8+djWZJS8P/vxn/+U0frz/smpjFMAiLWXBAjjnHNh7bz+XjQ0d6rMinnzSw7iNUQCLtISyMjj2WJ/9cPvtkJ8fu6LMdfbZvhP02Wf7xqRtiAJYpCVcdhm88oof9aahh83Ly4N77oGuXeE732lT2xcpgEWa24MP+lq/p53mXbDUr18/uPdeeO89OPnkNjM1TQEs0pxefRX+539gr718zq803De+4WtFPPigf4JoAwpiFyCSM95/33cBHjQIHn0U2rePXVH2ueAC+OADXyti8GA4/fTYFbUoBbBIc3jvPe/gzPwb/V69YleUnczgxht9BskZZ/iXl+PHx66qxWgIQqSpZszwb/ErK+HFF/WlW1MVFMCkSb5+8Kmn+gEsOUoBLNIUTz3l471m8MILMGxY7IpyQ7t2vnLc2LE+Pe3886GiInZVzU4BLNIYFRV+dNvhh/vBBK+/Dl/5SuyqcksawuedB9dfD4ce6ts55RAFsMiWmjfPhxwuu8znrb70EgwcGLuq3JSf77NJbr3V3+cRI+CRR2JX1WwUwCINtWaNr2o2bJivY3vnnXDffdC5c+zKct8pp8C0af6L7r//28eHZ8+OXVWTKYBF6rNihU+L2nZbX8f229+GmTN9vq9WN2s9w4bBlClwzTXw8st++cQT/d8iSymARepSXg7//CecdJIfpXXJJbD77n548b33asghlqIinys8ezace66PEe+8M+y3n293X1ISu8ItYqENLgEnGxs1alSYOnVq7DLiqajwAyleecVnMzzzDHzxhQ8vfO97Ph1q5MjYVUptS5f6KmoTJ8KcOT5mvN9+MGaMn48cCZ06tdjLm9m0EMKoRj9eASzQBgK4tNS7o6VLYfFin+g/fz7MnevBO2MGrF3r9+3bFw480Mcav/WtFv0BlmYSArzxBjz8MPz97779E/gQ0Q47+AyV7beH4mI/UrFfP98YtWdP//dt5FCSAliaxajevcPUsWNb/oU29f8tvb7m7SFsfKqs9POKiupTebmfyspg/Xo/rVvnp9WrYdUq/3NdBg6EHXeE4cO9Wxo92n9gNbab3ZYt83U5pk/3L0zff99/2ZaWbnzfggLo0sVPHTtChw4+Ba5dOx/yKCz0++TnV5/694drr1UAS/MYVVQUprbWVumbCrf0+pq3m214ysvz85o/DPn5/kNSWFj9g9Ohg5+6dPFhhO7d/dS7N/Tp48E7cKDfV9qGykr/9PPxx36+ZIl/Kiop8V/Sq1b5p6B166p/kZeW+i/28vLqX/iVlb679TPPNDmAtRaEuF12gVweghDJy/Ohh379YldSRbMgREQiUQCLiESiMWABwMxWAe/HrqMevYDPYxdRj2yoEbKjzmyocccQQpfGPlhjwJJ6vylfJrQGM5uqGptHNtSZLTU25fEaghARiUQBLCISiQJYUrfGLqABVGPzyYY6c75GfQknIhKJOmARkUgUwCIikSiA2zgzO9jM3jezuWZ2cex6UmY2yMxeMLOZZvaumZ2bXN/TzJ41sznJeY8MqDXfzKab2RPJ5W3NbErynt5vZkWR6+tuZpPM7D0zm2Vme2ba+2hm5yf/zjPM7D4za58J76OZ/dnMlpjZjBrX1fnemftdUu/bZlbv+qUK4DbMzPKBm4BDgGHAcWaWKdv6lgMXhBCGAXsAZya1XQw8H0LYHng+uRzbucCsGpevAq4LIQwFSoDvR6mq2g3A0yGEnYBd8Voz5n00swHAOcCoEMJwIB84lsx4H/8CHFzruk29d4cA2yen8cAt9T57CEGnNnoC9gSeqXH5EuCS2HVtotZHgYPwo/X6Jdf1ww8giVnXwOSH8BvAE4DhR28V1PUeR6ivG/AhyRfuNa7PmPcRGAB8AvTEDw57AvhWpryPQDEwo773DvgjcFxd99vUSR1w25b+x08tSK7LKGZWDHwVmAL0CSF8lty0COgTqazU9cCFQGVyeStgeQihPLkc+z3dFlgK3J4Mk9xmZp3IoPcxhLAQuAb4GPgMWAFMI7Pex5o29d5t8c+TAlgympl1Bh4CzgshrKx5W/A2I9o8SjM7HFgSQpgWq4YGKABGAreEEL4KrKHWcEMGvI89gLH4L4v+QCc2/tifkZr63imA27aFwKAalwcm12UEMyvEw/eeEMLfkqsXm1m/5PZ+wJJY9QF7A//PzD4C/ooPQ9wAdDezdJ2V2O/pAmBBCGFKcnkSHsiZ9D4eCHwYQlgaQigD/oa/t5n0Pta0qfdui3+eFMBt2xvA9sm3zUX4Fx+PRa4J8G+UgYnArBDCtTVuegwYl/x5HD42HEUI4ZIQwsAQQjH+3v0zhHA88AJwdHK32DUuAj4xsx2Tq8YAM8mg9xEfetjDzDom/+5pjRnzPtayqffuMeDEZDbEHsCKGkMVdYs18K5TZpyAQ4HZwAfAZbHrqVHXPvhHu7eBt5LTofgY6/PAHOA5oGfsWpN69weeSP48BHgdmAs8CLSLXNtuwNTkvXwE6JFp7yPwC+A9YAZwF9AuE95H4D58XLoM/zTx/U29d/gXsDclP0vv4LM6Nvv8TToU2cwOxj9y5QO3hRB+3egnExFpYxodwMkc0tn41KAF+MfZ40IIM5uvPBGR3NWUBdm/DswNIcwDMLO/4t9kbjKAe/XqFYqLi5vwktJUy5f75q+DBm14/bRp0z4PIfROLx+U9x2t0iRSw7OVD25iO+/Ga0oA1zXnbXTtO5nZePyoEAYPHsxU7bwb1YUXwo03brwBspnNj1ORSNvV4rMgQgi3hhBGhRBG9e7du/4HSIsqK4PCwthViAg0LYAzeg6p1K20FIqiLg0jIqmmBHDGziGVTVu7Fjp1il2FiEATxoBDCOVmdhbwDD4N7c8hhHebrTJpEStXQufOsasQEWjitvQhhCeBJ5upFmkFy5bBVlvFrkJEQIcitzkLF0LfvrGrEBFQALcpFRXw0Uew3XaxK4nAbMtOIq1AAdyGzJgB5eUwfHjsSkQEmjgGLNnlX//y8/32i1tHRrB6eo+0CQ6VG17fhLVTRGpTB9yG3H03DBu28WHIIhKHOuA2YvJkeOMNPwy5TUo713R8N+1sk07Y8myDy6SXK2t1vMnjQnq9OmRpAnXAbUBpKZx+OvTpAyeeGLsaEUmpA85xIcBll8Fbb8Ejj0DXrrEriqx2J1xb0vlaenthfp33t+R5QkWFX1FZ67I6Y2kAdcA5LARf/eyaa7wDHjs2dkUiUpM64By1YgWccw7ceSeceSb87nexK8owVR1pMqabNK5Vkx/yvfNNO2FLLlN1nowdV2zY6YbSUv9D0gmH9PZ07DjtkDeqQ9oidcA56JFHfLbD3XfDz37mX7zl6V9aJOOoA84hr70GV14Jjz8Ou+ziQfy1r8WuKsPV6kDTDrVqxDfpgEPehrMlrCD50WlfsMH9LC9Z6Wi9d8KhrMwvl5X75fLkvOpy2WbrkdymvijLVVTAww/D3nvDXnvBv/8NV1zhO14ofEUymzrgLDVrFtx/P9xzD8ydC8XFPs77v/+r5SabpNbshnS+r9XuTPOSMeLCpENu76vch3Z+bp07bnj/Uu90875MOuNkrDisXuPntTrjqlkU6ohzmgI4i8ye7aH7wAO+roOZH1Z8+eXw7W9Dgf41RbKKfmQz2OrV8PLL8Nxz8Oyz8M47Hrr77ONfrB11FPTrF7vKHJV2nsn0iJAO1aazGKrm/9YaMy7wzriiS3u/uZ3/iIWCZIy43J83f7V3wHlruvj5+lod8bov/fmSseXKtWv9YjILI+2YJbspgDNIWRlMmeKB+/zzfvhwebnv4bb33nD99XD00TBgQOxKRaQ5KIAjWrLEAzc9vfYarFnjTc/uu8MFF8CBB3r4dugQu9o2Lul0Q0h63XS+b9XtG44V5yVH0KUdcHkHv1zWORk7rkjGjJNZE+1KvKMtWr7eH7/cO2FW+Xn6bXna+VqtyxupWvNCY8iZTAHcStav98OBJ0/2sJ08GT780G/Lz4cRI2DcOBgzBvbfH3r2jFquiLQCBXALWL0a3n7bA3f6dD9/++3qpmnAANhjDz88ePRo73a1U3GWqL0GRDqbIZm9kJdcn54Xrk87VB/rrUxmTazt5Z1weUe/vHqAXy5cU5Sce2fcYYkPPrdb4p1w/her/PXWrfPnS8eMa8+eSFd1C7WOvJOMogBuosWLNwza6dNhzpzqT349e8JXvwrnnuthO3o0DBwYt2YRyQwK4AZatw5mzvSZCDVPixZV36e42MP2+ONht938zwMHaouxnJR2wukYbPKPXLE6mRWxPhnL/dLPi5JZDnnru/nj8nye8LqkU10z0B+/ro8/XaU3whSsbgdAx0V+RefPfDm7dp/78xUuWuF3LPHzsMZnS6RrUAStypbRFMC1VFTAvHkbB+3cuVCZ/F9u397XWvjWt6qDdtddoXv3uLWLSHZp0wG8ZMnGQfvuu5BMucTMdxAeMQKOPdbPR4yAoUOrF8USATY5b7hylY/ZWjJWnL/Gx267rPBOuHCQn+eVe4e7clvvhMt6JZ11X++gVwzxi8tX+P2KlnkH3elTnx7T4/0efv2S1f56y5b7A5JOvDJ53aq1J9Ix4kqNEcfUJgJ47VoP1tphu2RJ9X169/ZwPeWU6qDdeWd9OSYiLSfnAviLL/yLsDffrD6fPbu6QWnf3rdlP+yw6qAdMcK36xFpNrXmDYf0SLZ0Hu+XfqRbh2SNiKISn/WQzn5YvbwQgFXbeafasb93tkMHfwpAaaV/BJv9qf/HXf4VHyvu/FEvALp+7B1xpw/9cfkl3olXLivx+tK1LtJV2cpqzWuWVpG1ARwCfPrpxmH78cfV9xk0CEaO3HD4YLvtNHwgIpkhawK4stKHEV56yddHeOkl+Oyz6tt32AH23NN3f/jqV/3Uq1e8eqWNqz31JZ01kYzJpjttVH62GIC85SsB6LbSO9eOyWyH9st8jPfz3Xwe8awy7x527OvjZ/8zfAoAPQt8PvBji3YBYO4HfQHoOsvHmLt+lHTYq3r78370hdexZJnXkc6WSNeaSOqUlpWxAVxWBtOmVYftK69ASfLpacAAP1psjz28w911V+jSJWq5IiJbLKMCOARfUPzOO33JxZXeFLDDDr7c4r77+vKLxcWaWysZrp75tpVph5nMRqgaG07OC1f62O1WJd4Rt1vpnezqfv6t8DsjBgOwfDvvkE/e5hUAbhr6VwA+KfYO+t5d9gTglY+3BSD/Le9UenTfGoDO870zzv/UO+GwxjvpymRti3SX6Koj7TRrolnVG8BmNgi4E+gDBODWEMINZtYTuB8oBj4CjgkhlDSmiHnz4K67PHjnzfOZB0cdBUcc4Usv9u3bmGcVEclsDemAy4ELQghvmlkXYJqZPQucBDwfQvi1mV0MXAxctKUF/Pa38KMfeUc7ZgxMmODdrqZ/SU6rPW843QAjOZLNkrUebG0yb3i5z2Lo1MM72y4L/HzJ7r4g9C8XHw7Af+/8FgD7dp0NwMTB/wbg7b7PAfDn7fYB4MnZOwPwxQzviDsu9k6413+SWRNLfB5xWOEfQ4Ml84nXa6eO5lRvAIcQPgM+S/68ysxmAQOAscD+yd3uAF5kCwP41ls9fI86Cq67zmctiIi0FRa24DeZmRUDLwHDgY9DCN2T6w0oSS9vyqhRo8LUqVMBX4pxu+18JbBXXvFFxyUeM5sWQhiVXj4o7ztqcTJButtyMjshv4+P3ZLs0lw22Kf6rOvj84C/2Mnvt24nn2d88m6vAnBA55kADCv068vwf94Ji8YAMPmzbQBY/9pWAPR8zzvzDov8/oWLk7UmkjUn0nnE6U4dVavD5XBn/Gzlg83+zVODd0U2s87AQ8B5IYSVNW8LnuJ1vvNmNt7MpprZ1KVLl1ZdP3AgHHCAz9998snGFS8iks0a1AGbWSHwBPBMCOHa5Lr3gf1DCJ+ZWT/gxRDCjpt7npodMMCqVXDQQfDGG76wzbhxMHasH60mrUsdcIZLOuG8dt7ppvN1rSjdhdm/NKno67MmVg/2tSLSjthGeud6+JB3ATik69sA7N/Bx3T/sdaPvHttzfYA3D9npL/MFB9r3mqmH7HXYUGy1sSiZNZEMkYdklXfqtaayMFOOEoHnAwvTARmpeGbeAwYl/x5HPDolr54ly7w9NNw0UW+NsOxx/qMh/HjfViisrL+5xARyVb1dsBmtg/wMvAOkEbipcAU4AFgMDAfn4b2xeaeq3YHXFNFBbz4ItxxBzz0kC+g06OHT0Pbbz+fAzxyJBQWbslfTxpKHXCWSCfA24a9U177DTtjhviq/+medKuGeIchLEvvAAAN7klEQVS8Ykgy73iUd8TD+/jhpIdu5R3x19r7sfxPrR7u54t9tsQH7/UHoNtMf/7uc33tiA4f+WwJW5GsNZGuR5zM4giVoXqKRypLu+OW6IAbMgvi39TYdbuWMc1VSH6+T0MbMwZuvhkefRReeMGPhHv8cb9Px45+uHF6QMbXv67paiKSvbZoFkRTba4D3pxFi/wIuXQdiP/8x3+JmsGOO/q6DyNHVq8BoQ0tt5w64CyzqV2P07Hijj4GbElnbMnlij4+UWntQL+8dDfvwb7s60e6FW/na1Mc1Oc9v1+yNcfSUp8n/MYiPwJvzX/8h6xrsrFs99ne8RYuS3ZzXvS5l7duXdWuISE9ui5Lj6aL0gFngr594eij/QSwfDm8+qp/eTd9uofzffdV33+bbaoDeeRI37Wif38dviwimSUrAri27t3h0EP9lPr8cw/j9PTmm/DII9UNQo8e1UtS7rKLnw8frkV8JEul/7Frd8LJeWWypkO6vYsluyfnJTtzdFngj+u40OcVf9nbpx6tHORjvbdv5+cF2yZrUnTxxw/p4bMfPh7p47pLe/msi/XdvaPu+nGyh11yuXDxCsIXJUlNybhwqHU0XV6t9WGztENujKwM4Lr06uVT2g46qPq6Vat8uOI///FZFm+/7etNJLvEAL6wT82F2UeM8MV/9GWfZIX6hhCTL+tCErwVSz1ArdB/9G1lcohz1+SQ5A/8vPscX/ynZEcfeijp4dd/2s+XsyzcOjlUupMPL6z6itdR3tF/cNZ3SQK5YwHt2vl1eV8kB3Os8lCvWgw+nbpWVfMmhldyUM4EcF26dPFZFPvsU31dCDB//sbbEz31FKQb3BYVwU47bRzM2uFYRJpTTgdwXcy86y0u9tXWUuvXw/vve5echvK//gX33FN9n+7dfdii9lBG166t/bcQaaDaH+fTxX/WJ+e1O4rl3qW2X+8dc99PPCLK+ntHvL67d7Nr+vj0o7X9k01Eu/iwQnln71rX9k0WnC8sonMHH6bo0NmHOfI/9yU0Q4lPYatcV2vJy3SIog10wm0ugDelXTsP1F122fD6khKYMWPDbvmee6rXKgbfJTndnj4979tX3bKIbJ4CuB49evi84333rb4uBPjkE++W33rLT2++CZMmVd9n662rA3m33XzRoaFDFcqSYdIv7dIv6wp92ln5xwuTyx4RRav99qICv9xpK//Yt76Xd7Ol3f36dT39P3i5DwFT0d5Y393HoSsL/cr2Hf2+he39tfKTpTYrk/Hoqmlrpbm/UagCuBHMYPBgPx1+ePX1K1Z4KE+f7qE8fTpce61vrwQ+P3n0aD/tsYcfSNKjR5y/g4jEpwBuRt26bdwtl5bCzJkwdSpMmQKTJ/v6F+mw1g47eBiPHg177eVDIHkNXqNOpHnV3p4+JM1D+WJfyTCvKJkelIwVd1jsY8HtO3kn3Lmbd7mlPXy8t7xzPmUdkpkYyWyzsk5J7GydbBSa/IdPn7t6EfhkHLk0XeAn9xaDVwC3sKKi6mGIH/zAr1u50gN58mQP5aef9ulx4NPpDjgADjzQD8seMkTDFiK5SgEcQdeu8I1v+Amqp8a99BI8/7yfHnzQbysurl4j4+CDNWQhrazWLIrKL5PLSVeaLrpjK3w8N6/EO98Oi5OZDh3bU9HTu+TKwmTL+4Kko6hIxp87Jo9NXsPSQ5aThYXy8pOlLtPXrFr8Pfs7YgVwBqg5Ne7EE/3/0/vvV4fxQw/BxIl+cMg3vwnHHOPrJnfrFrtyEWkKBXAGMvMDQXbaCc4805fqfOMND+IHHoC//92HNg4+2MP4qKO0iL20sspam4mmXWmyMHtesiwl+fkUrPSxXpIZFKFDsnRmOuabv+F835DMjqgaeUvH4NIvR5LZEVXzhsneTlhf92SB/Hz/ou43v4GPPoLXXvNgnjYNTjjBFx/61a9g2bLYlYrIlsiK5SilbpWVvmbyb3/rh1J36AAnnww//rGH8pbQcpTSLGp9Y2z5+dWLx+clS2V2SD6updsqpQuvJB3yRkfApfOCk/mc6VhwOr+zat5wC28MGnVTTsk8eXn+5dyTT1Zv6XTrrTBsGFx/fdXGtSKSoRTAOWL4cPjzn2HOHNh/fzj/fJ9X/O67sSuTNiWEDU6hvJxQXuYrnlVUQEUFFStXU7FyNZXpafkKP5Usp7JkOWHVKj+tWeun8nLvcisDVAbMDDPzDrrGyZITZlkzd1MBnGO22QaeeALuvRfmzfMQfuWV2FWJSF0UwDnIDI47zten6NvXp649/3zsqqTN2qAbLvcZFJUVhLJSQlkplaVlflr3pZ/WrKNyzTrfzmjdOsLqNX5av95PFRU+3ps8b1VHbHlgeRt3whncESuAc9igQX5wx5Ahvp3TwoWxKxKRmhTAOa5PH3j4YV/v+JRTsnKqpOS6pCOu6oyTMePK9ev9lHTIobQ0OZX5qazcTyFQ52yupCOuvpx5nbACuA0YOhSuuMKnqr3wQuxqRCSlAG4jTjvNl8O8+ebYlYjUo9ZMiqrOOBn7TTvkqlPaCae3VwZCZXZ81FMAtxHt28O4cfDoo5BumCsicSmA25BvftMPKpo8OXYlIo1QuzOu1SETKus+1ZZBY8EK4DZkzz39/PXX49YhIk6robUh3br5rIgPPohdiUgLyMIpPuqA25jiYvj449hViAgogNucXr20bKVIplAAtzE9ekBJSewqRAS2IIDNLN/MppvZE8nlbc1sipnNNbP7zayo5cqU5tKxIyTbeIlIZFvSAZ8LzKpx+SrguhDCUKAE+H5zFiYto317BbBIpmhQAJvZQOAw4LbksgHfACYld7kDOLIlCpTmVVhYtZGAiETW0A74euBCqna/YytgeQgh3RVvATCgrgea2Xgzm2pmU5cuXdqkYqXpCgqqdngRkcjqDWAzOxxYEkKY1pgXCCHcGkIYFUIY1bt378Y8hTSjvDzfS05E4mvIgRh7A//PzA4F2gNdgRuA7mZWkHTBAwGtNpsFFMAimaPeDjiEcEkIYWAIoRg4FvhnCOF44AXg6ORu44BHW6xKaTZmWXnAkEhOaso84IuAH5rZXHxMeGLzlCQtSQEskjm2aC2IEMKLwIvJn+cBX2/+kqQlZcgiUCKCjoQTEYlGAdzGqAMWyRwKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJpEEBbGbdzWySmb1nZrPMbE8z62lmz5rZnOS8R0sXKyKSSxraAd8APB1C2AnYFZgFXAw8H0LYHng+uSwiIg1UbwCbWTdgP2AiQAihNISwHBgL3JHc7Q7gyJYqUkQkFzWkA94WWArcbmbTzew2M+sE9AkhfJbcZxHQp64Hm9l4M5tqZlOXLl3aPFWLiOSAhgRwATASuCWE8FVgDbWGG0IIAQh1PTiEcGsIYVQIYVTv3r2bWq+ISM5oSAAvABaEEKYklyfhgbzYzPoBJOdLWqZEEZHcVG8AhxAWAZ+Y2Y7JVWOAmcBjwLjkunHAoy1SoYhIjipo4P3OBu4xsyJgHvC/eHg/YGbfB+YDx7RMiSIiualBARxCeAsYVcdNY5q3HBGRtkNHwomIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIGhTAZna+mb1rZjPM7D4za29m25rZFDOba2b3m1lRSxcrIpJL6g1gMxsAnAOMCiEMB/KBY4GrgOtCCEOBEuD7LVmoiEiuaegQRAHQwcwKgI7AZ8A3gEnJ7XcARzZ/eSIiuaveAA4hLASuAT7Gg3cFMA1YHkIoT+62ABhQ1+PNbLyZTTWzqUuXLm2eqkVEckBDhiB6AGOBbYH+QCfg4Ia+QAjh1hDCqBDCqN69eze6UBGRXNOQIYgDgQ9DCEtDCGXA34C9ge7JkATAQGBhC9UoIpKTGhLAHwN7mFlHMzNgDDATeAE4OrnPOODRlilRRCQ3NWQMeAr+ZdubwDvJY24FLgJ+aGZzga2AiS1Yp4hIzimo/y4QQvg58PNaV88Dvt7sFYmItBE6Ek5EJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCuI3ZYw8499zYVYgIgIUQWu/FzJYC81vtBWVLbBNC6B27CJG2pFUDWEREqmkIQkQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hE8v8BHZ1UcAkE23QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% Example with entropic regularization\n",
    "\n",
    "\n",
    "def f(G):\n",
    "    return np.sum(G * np.log(G))\n",
    "\n",
    "\n",
    "def df(G):\n",
    "    return np.log(G) + 1.\n",
    "\n",
    "\n",
    "reg = 1e-3\n",
    "\n",
    "Ge = ot.optim.cg(a, b, M, reg, f, df, verbose=True)\n",
    "\n",
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, Ge, 'OT matrix Entrop. reg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD with Frobenius norm + entropic regularization\n",
    "-------------------------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "    0|1.693084e-01|0.000000e+00\n",
      "    1|1.610121e-01|-5.152589e-02\n",
      "    2|1.609378e-01|-4.622297e-04\n",
      "    3|1.609284e-01|-5.830043e-05\n",
      "    4|1.609284e-01|-1.111407e-12\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XeYVdW5x/HvO42hF0Gkj1ijYiwkYL0qds3VG0s0JmJi7LFFY0vTGI0aY4lRE68m0Vhii5oQy1VjSVRQEKMoKoggoAhIr9PW/eNdhxmGgZlhyjrnzO/zPPs5c/bZs/c7e+a88561117LQgiIiEjbK0gdgIhIe6UELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwSAPM7Hdm9pPUcbQmM1tmZkMTx/AnM/tFyhjamhKwbBQzO8nM3jGzFWY2x8xuN7Me8bXfxTf0MjMrN7OKWs+fauE49jGzWS25z7pCCKeHEK5szWO0FjN70cy+19B2IYQuIYRpG7H/P8Xf8bJayzc2Ltr2RwlYmszMLgCuBX4IdAdGAkOAZ82sJCasLiGELsDVwIOZ5yGEQxLEW9TWx6wnhsvN7PLUcdTVQufmulq/3y4hhAdb6TiNlg2/88ZQApYmMbNuwBXA2SGEp0MIFSGE6cCxQBnwrY3c7+Fm9paZLTKzV81sx1qvTTezC83sbTNbbGYPmlmpmXUGngL616q++sdk94iZ3WtmS4CTzKyDmd1kZp/G5SYz6xD3v4+ZzTKzy8xsfjzeCbWOv9ZHYzM7Isa6xMw+MrODN+Zn3sC5CGZ2pplNMbOlZnalmW0Rz8sSM3vIzEritj3NbIyZzTOzhfHrgfG1q4C9gN/Gc/PbWvs/y8ymAFNqrdvSzEriz3Z2XF9oZq+Y2U834ueYbmYXm9nbwHIzKzKzL8WqfJGZvWtm/13n23qb2bPx537JzIY08lgnxThvNLMvgMvj+u+a2eR4bp6pvT8zO9DMPoh/U7fF4zX4aaFFhRC0aGn0AhwMVAJF9bx2N/BAnXWXA/c2sM+dgbnACKAQGA1MBzrE16cDrwP9gV7AZOD0+No+wKx6jlkBHIkXGR2BnwNjgU2BPsCrwJW19lEJ3AB0AP4LWA5sE1//E/CL+PVXgcXAAXHfA4BtG3HeLgcub+Q5DsATQDdge2A18DwwFP/E8R4wOm67CXAU0AnoCjwMPF5rXy8C36tn/8/Gc9mx1rot49c7AAuBLwE/iuetcD2xrjk39bw2HXgLGBR/B8XAVOAyoATYD1ha5zwvBfaOv4ebgX838pydFH+HZwNF8XhHxON9Ka77MfBq3L43sAT4enzt3Pg3873GHK+lFlXA0lS9gfkhhMp6Xvssvt5UpwK/DyGMCyFUhRDuxpPOyFrb/CaE8GkIYQHwd2CnBvb5Wgjh8RBCdQhhJXAC8PMQwtwQwjy8iv92ne/5SQhhdQjhJeAfeFVf18nAH0IIz8Z9zw4hvN/0H7lB14UQloQQ3gUmAf8XQpgWQliMV/07A4QQvgghPBpCWBFCWApchf8DacgvQwgL4rlZSwhhEvAL4HHgQuDbIYSqDezrwljRLjKz+XVe+00IYWY8zkigC3BNCKE8hPBPYAxwfK3t/xFCeDmEsBpP/ruZ2aBG/DwAn4YQbgkhVMbjnR5/zsnx7/VqYKdYBR8KvBtC+Gt87TfAnEYep8UoAUtTzcc/JtbXxtYvvt5UQ4ALar2JF+FVU/9a29R+c6zA38gbMrPO8/7AjFrPZ9TZ/8IQwvINvJ4xCPiogWMDEJsDMj/PJcAltX7GMQ18++e1vl5Zz/Mu8RidzOz3ZjYjNre8DPQws8IG9l/3/NR1N/57eTKEMKWBba8PIfSIS91/wLWP0x+YGUKorrVuBv4pYp3tQwjLgAXU/3uoT92faQhwc63fwQLA4vH61zlWAFr1Ym59lIClqV7Dq9Ov115pZl2AQ/CPyk01E7iq1pu4RwihUwjhgUZ87/qG86u7/lP8DZkxOK7L6BnblNf3eu1Yt2hEXIQQDs/8PMA1eOWX+fkOb8w+GuECYBtgRAihG/7xHTzRQOPPT1234dXpQWa2ZzPiq32cT4FBZlY77wwGZtd6vqbajX9Tvaj/99DQscB/V6fV+bvqGEJ4Ff+0NrDWsaz287aiBCxNEj8CXwHcYmYHm1mxmZUBD+EVxJ83Yrf/C5xuZiPMdTazw8ysayO+93NgEzPr3sB2DwA/NrM+ZtYb+Clwb51trogXofYCDsfbU+u6C/iOmY0yswIzG2Bm2zYiztbSFa+IF5lZL+BndV7/HG87bjQz+zawK96ueg5wd0yGzTUO//RyUfy72Qf4GvCXWtscamZ7xouMVwJjQwgNVevr8zvgUjPbHsDMupvZMfG1fwDDzOzI+GnuLGCzjTzORlMCliYLIVyHX0i5Hr+QMQ6vNkbFtrum7m88cArwW/ziz1T8zd+Y730fT67T4kfN9X1c/QUwHngbeAd4M67LmBOP/SlwH36Rb5223RDC68B3gBvxi3EvsXZl3dZuwi84zccvlj1d5/WbgaNjL4DfNLQzMxsc93liCGFZCOF+/Lzd2NxAQwjleMI9JMZ7WzxO7fN8P/5PZAH+T2BNr5rYa+IEGimE8BjeXfIvsXlmUjw2IYT5wDHAdcAXwHb4z9nkv9/msHhFUKTdipXYvSGENv8IKtkhNovMAk4IIbzQVsdVBSwi7ZKZHWRmPcz7g1+Gt5uPbcsYlIBFpL3aDe/RMh9vGjmyvm55rUlNECIiiagCFhFJJCcGrJDW17t371BWVpY6DJGcMmHChPkhhD4b+/1KwAJAWVkZ48ePTx2GSE4xsxkNb7V+aoIQEUlECVikPVm8GF59FWbMAF2AT04JWCTfhQAPPww77AA9esAee0BZGWy6KVxyCSxdmjrCdksJWCSfzZkDo0bBscdCYSFcdRU88QTcdhvsuy9cey1ssw0880zqSNslXYQTyVezZnnynTULbr0VTjvNk3DGGWfAuHFwyinwta/BAw/AUUeli7cdUgUsko9mz4a99vIK+P/+D848c+3kmzFiBPzrX/CVr3iV/HB9A8BJa1ECFsk3FRWeTOfPh+ef9zbfDene3ZsgdtsNRo+GSZPaJk5RAhbJOxdd5D0d7roLhg9v3Pd06QKPPOLJ+KijYMmS1o1RACVgkfwyZgzcdBOcc45XwU2x2Wbw4IPw0Udw1lmtE5+sRQlYJF8sW+ZtvTvsAL/61cbtY++94bLL4N574dlnWzY+WYcSsEi+uPxymDkTfv97KCnZ+P1cdhlstZUn81WrWiw8WZcSsEg+ePttb3o45RTYfffm7au01PsJT50Kv/xly8Qn9VICFskHF10E3brBNde0zP723x++8Q1vyvi0sZMSS1MpAYvkuhde8G5kl10GvXq13H6vvhoqK+HnP2+5fcpalIBFclkIPp7DwIHw/e+37L6HDvW75+68Ez78sGX3LYASsEhue+wxeP11uOIKb7ttaT/+se/3Jz9p+X2LErBIzgoBfvEL77Fw4omtc4y+feG88/wW5cmTW+cY7ZgSsEiueuopmDgRLr0UilpxXK3zzoOOHVvuAp+soQQskosy1e/gwfCtb7XusXr3hlNPhfvug2nTWvdY7YwSsEgueukleO01735WXNz6x7vwQh9N7brrWv9Y7YgSsEguuu46n9Hiu99tm+MNGOAjpf3pTzB3btscsx1QAhbJNe+95+2/3/++t822lR/8AFav9rvkpEUoAYvkmptu8q5hp5/etsfddls47DBPwCtXtu2x85QSsEgumTsX7rnHu5316dP2x7/gApg3zy/ISbMpAYvkkt/9zpsBzj8/zfH32Qd22gluuEHT2rcAJWCRXFFeDrffDocc4s0BKZh58p882ac7kmZRAhbJFY884pNsnn122ji+8Q1v/rjllrRx5AElYJFcccstftvxQQeljaNDB78x4+9/h48/ThtLjlMCFskF48fD2LHe9awgC962Z5zhcahLWrNkwW9SRBp0yy0+c/FJJ6WOxA0Y4LMn33knrFiROpqcpQQsku3mz/fZik880We9yBbf/z4sWgQPPJA6kpylBCyS7e66y7uenXlm6kjWtueeMGwY3HqruqRtJCVgkWxWVeVdz/bZB7bfPnU0azPzfwoTJ3r7tDSZErBINnvySZgxA846K3Uk9fvWt7xZ5Le/TR1JTlICFslmt94K/fvDEUekjqR+Xbr4KGkPPwyff546mpyjBCySraZM8dmOTzutbcb83VhnngkVFd4jQppECVgkW91+u081dMopqSPZsG23hVGj4Pe/92nspdGUgEWy0YoV8Mc/el/bfv1SR9Ows86CmTNhzJjUkeQUJWCRbHT//d7HNlsvvtX1ta/BoEHeZi2NpgQskm1C8EQ2bJj3tc0FRUU+QPxzz8H776eOJmcoAYtkm1degbfe8jvNzFJH03jf+x6UlKhLWhMoAYtkm1tugR494IQTUkfSNJtuCscdB3ffDUuWpI4mJygBi2ST2bPh0Ufh5JOhc+fU0TTd2WfDsmU+e7I0SAlYJJv87ndQXZ194z401vDhMGKEN0NUV6eOJuspAYtki1WrvC/t4YfD0KGpo9l4Z5/tN5E8/XTqSLKeErBItrjvPp9x+LzzUkfSPMcc432Xb7wxdSRZTwlYJBuE4Alrxx1h331TR9M8JSXeg+O55+Cdd1JHk9WUgEWywXPPwbvv+ozDudT1bH1OOw06doSbb04dSVZTAhbJBjfeCH37wvHHp46kZWyyic/gce+9MHdu6miylhKwSGrvvANPPeW3HXfokDqalnP++VBerunrN0AJWCS1667zPr+5Mu5DY22zjY9jfOut3jdY1qEELJLSjBk+qeUpp0CvXqmjaXkXXwwLF8L//m/qSLKSErBISjfc4BfdfvCD1JG0jpEj4b/+y3/O8vLU0WQdJWCRVObO9crwhBN8KMd8dfHFMGsW/PnPqSPJOkrAIqn86lc+3fyll6aOpHUdfDDsuitcdZVPXSRrKAGLpDB3Ltx2G3zzm36xKp+ZweWXw8cfqwquQwlYJIXrr/exH37849SRtI3DDlMVXA8lYJG2NmeOd806/vj8r34zMlXwtGkaqrIWJWCRtnbFFd4j4PLLU0fStg47zHtFXH65TzoqSsAiberDD73nw2mnwZZbpo6mbZn5TSeffqoxIiIlYJG2dNllPkjNT3+aOpI09trLZ1C+5hr44ovU0SSnBCzSVl5+2acbuvBCnz+tvfrlL/3W5J/9LHUkySkBi7SFykofI3fIEPjhD1NHk9b22/uUS7ff7rM/t2NKwCJt4dZbfdSzG2+ETp1SR5PelVf6kJVnndWu545TAhZpbZ9+6m2+Bx4IRx6ZOprs0KMHXHstvPqqT2PfTikBi7SmELzHQ3m5zxScD7NdtJTRo/2i3Pnnw+zZqaNJQglYpDXdey+MGQNXXw1bbZU6muxSUAB/+IP/czr1VP9n1c4oAYu0llmz4JxzYI89/FHWteWW3iviySc9GbczSsAiraGiAo47zns//PGPUFiYOqLsdfbZPhP02Wf7xKTtiBKwSGv40Y/glVf8rjc1PWxYQQHcdx906wbHHNOupi9SAhZpaQ8/7GP9nn66V8HSsH794P774f334bvfbTdd05SARVrSq6/Ct78Nu+/ufX6l8fbbz8eKePhh/wTRDhSlDkAkb3zwgc8CPGgQPPEElJamjij3XHABfPSRjxUxeDCccUbqiFqVErBIS3j/fa/gzPyKfu/eqSPKTWZwyy3eg+TMM/3i5amnpo6q1agJQqS5Jk3yq/jV1fDii7ro1lxFRfDIIz5+8Gmn+Q0seUoJWKQ5nnrK23vN4IUXYLvtUkeUHzp08JHjjjjCu6edfz5UVaWOqsUpAYtsjKoqv7vt8MP9ZoLXX4cvfSl1VPklk4TPOw9uugkOPdSnc8ojSsAiTTVtmjc5/OhH3m/15Zdh4MDUUeWnwkLvTXLHHX6ehw2Dxx9PHVWLUQIWaazly31Us+2283Fs77kHHngAunRJHVn+O+UUmDDB/9H9z/94+/CHH6aOqtmUgEUasnixd4vafHMfx/brX4f33vP+vhrdrO1stx2MGwfXXw//+pc/P/FE/13kKCVgkfpUVsI//wknneR3aV16Key6q99efP/9anJIpaTE+wp/+CGce663EW+/Pey9t093v3Bh6gibxEI7HAJO1jV8+PAwfvz41GGkU1XlN1K88or3ZnjmGViwwJsXvvlN7w61yy6po5S65s3zUdTuugumTPE24733hlGj/HGXXaBz51Y7vJlNCCEM3+jvVwIWaAcJuLzcq6N58+Dzz72j/4wZMHWqJ95Jk2DFCt92s81g//29rfGgg1r1DSwtJAR44w147DH4xz98+ifwJqKtt/YeKlttBWVlfqdiv34+MWqvXv773cimJCVgaRHD+/QJ4484ovUPtL6/t8z62q+HsO5SXe2PVVU1S2WlLxUVsHq1LytX+rJsGSxd6l/XZ+BA2GYb2GEHr5ZGjPA3rNp2c9sXX/i4HBMn+gXTDz7wf7bl5etuW1QEXbv60qkTdOzoXeA6dPAmj+Ji36awsGbp3x9uuEEJWFrG8JKSML6tpkpfX3LLrK/9utnaS0GBP9Z+MxQW+pukuLjmjdOxoy9du3ozQo8evvTpA337euIdONC3lfahuto//XzyiT/OneufihYu9H/SS5f6p6CVK2v+kZeX+z/2ysqaf/jV1T679TPPNDsBaywIcTvuCPncBCFSUOBND/36pY5kDfWCEBFJRAlYRCQRtQELAGa2FPggdRwN6A3MTx1EA3IhRsiNOHMhxm1CCF039pvVBiwZHzTnYkJbMLPxirFl5EKcuRJjc75fTRAiIokoAYuIJKIELBl3pA6gERRjy8mFOPM+Rl2EExFJRBWwiEgiSsAiIokoAbdzZnawmX1gZlPN7JLU8WSY2SAze8HM3jOzd83s3Li+l5k9a2ZT4mPPLIi10MwmmtmY+HxzMxsXz+mDZlaSOL4eZvaImb1vZpPNbLdsO49mdn78PU8yswfMrDQbzqOZ/cHM5prZpFrr6j135n4T433bzBocv1QJuB0zs0LgVuAQYDvgeDPLlml9K4ELQgjbASOBs2JslwDPhxC2Ap6Pz1M7F5hc6/m1wI0hhC2BhcDJSaKqcTPwdAhhW+DLeKxZcx7NbABwDjA8hLADUAgcR3acxz8BB9dZt75zdwiwVVxOBW5vcO8hBC3tdAF2A56p9fxS4NLUca0n1ieAA/C79frFdf3wG0hSxjUwvgn3A8YAht+9VVTfOU4QX3fgY+IF91rrs+Y8AgOAmUAv/OawMcBB2XIegTJgUkPnDvg9cHx9261vUQXcvmX+8DNmxXVZxczKgJ2BcUDfEMJn8aU5QN9EYWXcBFwEVMfnmwCLQgiV8Xnqc7o5MA/4Y2wmudPMOpNF5zGEMBu4HvgE+AxYDEwgu85jbes7d01+PykBS1Yzsy7Ao8B5IYQltV8LXmYk60dpZocDc0MIE1LF0AhFwC7A7SGEnYHl1GluyILz2BM4Av9n0R/ozLof+7NSc8+dEnD7NhsYVOv5wLguK5hZMZ587wsh/DWu/tzM+sXX+wFzU8UH7AH8t5lNB/6CN0PcDPQws8w4K6nP6SxgVghhXHz+CJ6Qs+k87g98HEKYF0KoAP6Kn9tsOo+1re/cNfn9pATcvr0BbBWvNpfgFz7+ljgmwK8oA3cBk0MIN9R66W/A6Pj1aLxtOIkQwqUhhIEhhDL83P0zhHAC8AJwdNwsdYxzgJlmtk1cNQp4jyw6j3jTw0gz6xR/75kYs+Y81rG+c/c34MTYG2IksLhWU0X9UjW8a8mOBTgU+BD4CPhR6nhqxbUn/tHubeCtuByKt7E+D0wBngN6pY41xrsPMCZ+PRR4HZgKPAx0SBzbTsD4eC4fB3pm23kErgDeByYBfwY6ZMN5BB7A26Ur8E8TJ6/v3OEXYG+N76V38F4dG9x/s25FNrOD8Y9chcCdIYRrNnpnIiLtzEYn4NiH9EO8a9As/OPs8SGE91ouPBGR/NWcAdm/CkwNIUwDMLO/4Fcy15uAe/fuHcrKyppxSGmuRYt88tdBg9ZeP2HChPkhhD6Z5wcUHKNRmkRqebb64fVM573xmpOA6+vzNqLuRmZ2Kn5XCIMHD2a8Zt5N6qKL4JZb1p0A2cxmpIlIpP1q9V4QIYQ7QgjDQwjD+/Tp0/A3SKuqqIDi4tRRiAg0LwFndR9SqV95OZQkHRpGRDKak4Cztg+prN+KFdC5c+ooRASa0QYcQqg0s+8Dz+Dd0P4QQni3xSKTVrFkCXTpkjoKEYFmTksfQngSeLKFYpE28MUXsMkmqaMQEdCtyO3O7Nmw2WapoxARUAJuV6qqYPp02GKL1JEkYNa0RaQNKAG3I5MmQWUl7LBD6khEBJrZBiy55aWX/HHvvdPGkRWsoM7TtaveUF3nRsBQXee5bhSU5lMF3I7cey9st926tyGLSBqqgNuJsWPhjTf8NuR2KVOx1mnfXVP5Ziri+NwK63x/piKOlXCo83yd44g0girgdqC8HM44A/r2hRNPTB2NiGSoAs5zIcCPfgRvvQWPPw7duqWOKLE1FWqmkvUaxAoylayXvlYYS+BMxVxQp824Om5fVeX7qYr7i8/XVMaqiGUDVAHnsRB89LPrr/cK+IgjUkckIrWpAs5TixfDOefAPffAWWfBb36TOqIss75KmFjBZirfOHScFcW3SmGsWer2oqiO31fhs6iHyviYea7KWOqhCjgPPf6493a491746U/9wluBftMiWUcVcB557TX45S/h73+HHXf0RPyVr6SOKstlKtEQ23JDnV4Sa76IX8VK2ErioMrFRWu/HvdnlXF/FRW+fvVqf15eEderMhZVwDmvqgoeewz22AN23x3+/W+4+mqf8ULJVyS7qQLOUZMnw4MPwn33wdSpUFbm7bzf+Y6Gm2yWWIGuqUzXrPb1mYolZPoLxwo4lPoo96E49p7ItPlkKuJyr3gLVpX7+pWr/OVVmcdYIVfWrYxVEeczJeAc8uGHnnQfesjHdTDz24qvugq+/vU1n45FJEfoLZvFli2Df/0LnnsOnn0W3nnHk+6ee/qFtaOOgn79UkeZpzKVcKxIM3fCVcdHq1OZWqx4Q4m/pao6eRtxdQeviNdUzPH7C1Z7hVu4zCvfgqUrfbtlK/yxbmWcqYir167MJbcpAWeRigoYN84T7vPP++3DlZU+h9see8BNN8HRR8OAAakjFZGWoASc0Ny5nnAzy2uvwfLlXuXuuitccAHsv78n344dU0fbzlVneknE3gor174TLvNYkKmcY0WcqYAruvhjZWmslDNNxZX+iy1a4bcoliz1XhJFC70iLly83LdfsswfV/r66tibQr0ncpsScBtZvdpvBx471pPt2LHw8cf+WmEhDBsGo0fDqFGwzz7Qq1fScEWkDSgBt4Jly+Dttz3hTpzoj2+/7YPigDchjBzptwePGOHVrmYqzhF12oYzo6JZpgKO6wtjP1+r6urbmf+Cqzp4BVze2duEKzr5Y4h32BWW+1uyeFkpAKWLvDIuned/PEXzvRIuXLTEv29ZrJBX12krVkWcE5SAm+nzz9dOtBMnwpQpNX//vXrBzjvDued6sh0xAgYOTBuziGQHJeBGWrkS3nvPeyLUXubMqdmmrMyT7QknwE47+dcDB2qKsbyWaRsur45P/bnFNlqLd8J1WB37AceKGDoAUNXB/zhWdfe1Kztn2pD9eeFqbywuWdwJgNL53mbc+fOe/vwz7zVROH+xf99ir4yrM/2MKzNtxaqIs5EScB1VVTBt2rqJdupUyIxAWFrqYy0cdFBNov3yl6FHj7Sxi0huadcJeO7cdRPtu+/CCi8qMPMZhIcNg+OO88dhw2DLLf3CmcgadduGM5VwbAsuiBVpSezXW7jC/1sXlsfuLcHfiiv8hjrKe8e23P5ewa403/+y5d6/eNEXvn2nz7107vypV9adZ/txiud4RcxCf6zOtBVXZPo1qz9xNmgXCXjFCk+sdZPt3Lk12/Tp48n1lFNqEu322+vimIi0nrxLwAsW+IWwN9+sefzww5omsNJSn5b9sMNqEu2wYT5dj0iLyVTEFd57oSrTSyJ2hSlY7h+zuizxttzipV7JFsUSeFmFvzVXFfp+em3mlWzfft4Lgi394bMl3kvis8+9Ai6d7W3FXWZ6Zd31k15xvbcNM3+Rx7V0KVCrP7Eq4iRyNgGHAJ9+um6y/eSTmm0GDYJddlm7+WCLLdR8ICLZIWcScHW1NyO8/LKPj/Dyy/DZZzWvb7017Labz/6w886+9O6dLl6RtcQKs3q1X8m12FZssW24NI4BUbwo9m5Y5EPaLVnibb4LV/r66iHePWL4ZjMB2LfPBwCs3ty3e2ux93F8e3Z/ABbN8Da0LjP6ANB9urc9d/zEK+DCeQsACEu9sq6ObdSqiNtG1ibgigqYMKEm2b7yCixc6K8NGOB3i40c6RXul78MXbtucHciIlknqxJwCD6g+D33+JCLS2Kz1dZb+3CLe+3lwy+WlalvreSo9fSWyMyUURAr0W4LYiU8Lz7G/r9LvvA23ee38LbehZv7+sP7vA3AgQPfAWBBP6+gn996OwCem7kNADM/9gq4y8e+nx7TvA25U6yIiz73irh6SWwjVkXcqhpMwGY2CLgH6AsE4I4Qws1m1gt4ECgDpgPHhhAWbkwQ06bBn//siXfaNO95cNRR8LWv+dCLm222MXsVEcluFhq4Q8bM+gH9QghvmllXYAJwJHASsCCEcI2ZXQL0DCFcvKF9DR8+PIwfP36tdb/+NVx4oVe0o0bBiSd6tavuX23LzCaEEIZnnh9QcIxunUohfrSzEu8NUdAlvhE28Up41WCvYBcP9deXbOEvd9jKPy4eUvYeAMf1HAfANsXe5vxprLQfW7KTP878MgBzp24CQNdpfmW6x0demXea4fuzWBGvucOuHfeaeLb64Rb/3N1gBRxC+Az4LH691MwmAwOAI4B94mZ3Ay8CG0zAdd1xhyffo46CG2/0XgsiIu1Fk9qAzawM2BkYB/SNyRlgDt5E0WgffwxE55M6AAAPjklEQVSnnw7Dh8P99/ug4yLtXqaNOI5uVpW5ky72Gy5d6JVo6Wfehtt1trfhLprl/Yj/usUIAF7ayjsKHzX4LQCO7TYRgIs3mbLW84cG7QzAI0P9cVZZ3O9HsY34I7+6namIC+dmek2oH3FLaPSsyGbWBXgUOC+EsKT2a8HbMer9yGpmp5rZeDMbP2/evDXrBw6Efff1/rtPPrlxwYuI5LIG24ABzKwYGAM8E0K4Ia77ANgnhPBZbCd+MYSwzYb2U7cNeOlSOOAAeOMNH9hm9Gg44gi/W03altqAc0SBt9UWlHi/X+vuFXDo6xXriiGxIt7CP9wuHeqV6aZbfgHA/wz6jz9288q4f7wr6YMKr8X+stAr6Keme++J1VN8f90+8sN3n+Z38pV+4nfU8YVfd18z1kRm0Os8HH2tNdqAG6yAzcyAu4DJmeQb/Q0YHb8eDTzR1IN37QpPPw0XX+xjMxx3nPd4OPVU7/ebGX1MRCQfNaYXxJ7Av4B3gExKvAxvB34IGAzMwLuhLdjQvurrBZFRVQUvvgh33w2PPuoD6PTs6d3Q9t7b+wDvsgsUFzflx5PGUgWcozIVcamPL1zQzdtsqzMV8WB/vmioV8TLNve3cPehXrmOGvihP3bz3hO9Cr0f8nurfebXMfN2BODNjwcD0OEj/3jabZr/eXT72OeoK5nt+wuZirj2eMR5Ug2n6gXxb2B9Bx7VUoEUFno3tFGj4Lbb4Ikn4IUX/E64v//dt+nUyW83ztyQ8dWvqruaiOSuRrUBt5QNVcAbMmeO3yGXGQfiP//xf6pmsM02Pu7DLrvUjAGhCS2bThVwnqhTEVtXvyMu9PE3xcohXhEvHuIfJZcNiXPaDfE23B0HfArATt1nAdChwHs5fLjc74YaP8f7ii6e4b0uun7sx+s23duaO0+Pc9bNjZXwosV5M19dkgo4G2y2GRx9tC8AixbBq6/6xbuJEz05P/BAzfZDhtQk5F128Vkr+vfX7csikl1yIgHX1aMHHHqoLxnz53syzixvvgmPP17zz7Znz5ohKXfc0R932EGD+EieyYy6lpnWJY7lUBDvZOsUxw3uNM3vqFs1wHs5LB3sbXnvDtoagAkDNgegR1/v79uv25K1HgtiW/KCjl4Jl3fzinpVj9hrYoa3FXeY3QmL7cKZGZzVd7hGTibg+vTu7V3aDjigZt3Spd5c8Z//eC+Lt9/28SZiH3LAB/apPTD7sGE++I8u9kleyCTkVTHZZSYLjSNdlc71xFs63RNpj838cfkAT6DL+3vTxZS+fit05SZ+Y0hhZ99PYSd/vqqff7ysLvYmicqO3gTStVsvOs32AYMK53rXNctMk6SJQ/MnAdena1fvRbHnnjXrQoAZM9adnuippyAOUEVJCWy77bqJWTMci0hLyusEXB8zr3rLyny0tYzVq+GDD7xKziTll16C++6r2aZHD2+2qNuU0a1bW/8UIhspVsRhtT9WxRsnbIlfPCue6zds9JwRmxI29Yp4VV8f/nJ53zhVUm//iFjePVatHeLuY0ZZlZkMoaCIqhKvsjt18heLM13mFnglHJavfRNHrl+sa4p2l4DXp0MHT6g77rj2+oULYdKktavl++6rGasYfJbkzPT0mcfNNlO1LCIbpgTcgJ49vd/xXnvVrAsBZs70avmtt3x580145JGabTbdtCYh77QT7LqrJ2olZckqdScPje2xFi/i2QK/gNZ5tndn69TdL+JV9PHnq3v7KFqrenjbb0UX/wOvioNrVRfC6u6ZP3qvfEsL/Qbckg6efgoW+MZhWZwWKdM2HAciyueLdUrAG8EMBg/25fDDa9YvXuxJeeJET8oTJ8INN/j0SuD9k0eM8GXkSL+RpGfPND+DiKSnBNyCundft1ouL4f33oPx42HcOBg71se/yDRvbb21J+MRI2D33b0JpKDRY9SJtLD1TJmUqUoLFnm7bfFc79lQ0s0r4c7dvZ23sqevL+/qbcSVnQrWtAuHWAhXdPVq2YK3KxfHP/iC4lgRF3klHFb4bc753DasBNzKSkpqmiG+9z1ft2SJJ+SxYz0pP/20d48D7063776w//5+W/bQoWq2EMlXSsAJdOsG++3nC9R0jXv5ZXj+eV8efthfKyurGSPj4IPVZCFtLFNthvX0J443V9gCb98t/jxWxp29uq3uUkpVF3+tuoNXvqFg7YqiqlPsdB8r4oJMxRGHyiRWwpmbStZUwnnQNqwEnAVqd4078UT/m//gg5pk/OijcNddfnPIgQfCscf6uMndu6eOXESaQwk4C5n5jSDbbgtnneVDdb7xhifihx6Cf/zDmzYOPtiT8VFHaRB7aWOZ/sSZx0zviZVerVq89dlKO1Cc+ePs6I+hNPZ6KInpJ1PxZtqfO8TB5qvX/qO2uF11ZnCfeHGbUL3W9+cSXe7JAYWFfqHuV7+C6dPhtdc8MU+YAN/6lg8+dOWV8MUXqSMVkabIieEopX7V1T5m8q9/7bdSd+wI3/0u/PCHnpSbQsNRSovKVLVWgMW2XIu9HCxOp0SHOGRmZuCVosK1vzeTm2J/4JAZxCdTAWd6R9TpsdFalXCSKYkkexUU+MW5J5+smdLpjjtgu+3gppu86UJEspcScJ7YYQf4wx9gyhTYZx84/3zvV/zuu6kjk3YpBF+qqwgV5YSKcqpXrqR65UqqliyjaskyqhcsistCqhcsJCxY5MuSpb4sX+lLeYVXv6Hal8JCX4qLobgYKyrypbDQq+2CuJhlfR9OJeA8M2QIjBkD998P06Z5En7lldRRiUh91AsiD5nB8cf7MJz77+9d1/72N2+uEEmmTp/idXpQZPr9WkF8HuvDzPq6t4hmpk0vWHs7I7P/zPbZ20tCFXAeGzTIb+4YOtSnc5o9O3VEIlKbEnCe69sXHnvMLxyfckpWFgHS3sX24lBZGZcKQmUF1atX+7JyFdUrVxFWrvRl1WpfMttXVdX0gACvoK0AKzCswNY8z8Y2YSXgdmDLLeHqq72r2gsvpI5GRDKUgNuJ00/34TBvuy11JCINyPSgqNWTguqqNZVupkIOFZVrLVRV+RJ7S4TqQKiu5yNfFlXCSsDtRGkpjB4NTzwBcQYYEUlMCbgdOfBAn3h07NjUkYhshPVUxutUyJk24Uy/4cyShZSA25HddvPH119PG4eIOPUDbke6d/deER99lDoSkVaQg118VAG3M2Vl8MknqaMQEVACbnd699awlSLZQgm4nenZExYuTB2FiEATErCZFZrZRDMbE59vbmbjzGyqmT1oZiWtF6a0lE6dIE5aICKJNaUCPheYXOv5tcCNIYQtgYXAyS0ZmLSO0lIlYJFs0agEbGYDgcOAO+NzA/YDHomb3A0c2RoBSssqLoaKioa3E5HW19gK+CbgItaM68YmwKIQQmV8PgsYUN83mtmpZjbezMbPmzevWcFK8xUV+c0YIpJegwnYzA4H5oYQJmzMAUIId4QQhocQhvfp02djdiEtqKCgZhhVEUmrMTdi7AH8t5kdCpQC3YCbgR5mVhSr4IGARpvNAUrAItmjwQo4hHBpCGFgCKEMOA74ZwjhBOAF4Oi42WjgiVaLUlqMWU7eMCSSl5rTD/hi4AdmNhVvE76rZUKS1qQELJI9mjQWRAjhReDF+PU04KstH5K0piwZBlVE0J1wIiLJKAG3M6qARbKHErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCJKwCIiiSgBi4gkogQsIpKIErCISCKNSsBm1sPMHjGz981sspntZma9zOxZM5sSH3u2drAiIvmksRXwzcDTIYRtgS8Dk4FLgOdDCFsBz8fnIiLSSA0mYDPrDuwN3AUQQigPISwCjgDujpvdDRzZWkGKiOSjxlTAmwPzgD+a2UQzu9PMOgN9QwifxW3mAH3r+2YzO9XMxpvZ+Hnz5rVM1CIieaAxCbgI2AW4PYSwM7CcOs0NIYQAhPq+OYRwRwhheAhheJ8+fZobr4hI3mhMAp4FzAohjIvPH8ET8udm1g8gPs5tnRBFRPJTgwk4hDAHmGlm28RVo4D3gL8Bo+O60cATrRKhiEieKmrkdmcD95lZCTAN+A6evB8ys5OBGcCxrROiiEh+alQCDiG8BQyv56VRLRuOiEj7oTvhREQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUmkUQnYzM43s3fNbJKZPWBmpWa2uZmNM7OpZvagmZW0drAiIvmkwQRsZgOAc4DhIYQdgELgOOBa4MYQwpbAQuDk1gxURCTfNLYJogjoaGZFQCfgM2A/4JH4+t3AkS0fnohI/mowAYcQZgPXA5/giXcxMAFYFEKojJvNAgbU9/1mdqqZjTez8fPmzWuZqEVE8kBjmiB6AkcAmwP9gc7AwY09QAjhjhDC8BDC8D59+mx0oCIi+aYxTRD7Ax+HEOaFECqAvwJ7AD1ikwTAQGB2K8UoIpKXGpOAPwFGmlknMzNgFPAe8AJwdNxmNPBE64QoIpKfGtMGPA6/2PYm8E78njuAi4EfmNlUYBPgrlaMU0Qk7xQ1vAmEEH4G/KzO6mnAV1s8IhGRdkJ3womIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAYuIJKIELCKSiBKwiEgiSsAiIokoAbczI0fCueemjkJEACyE0HYHM5sHzGizA0pTDAkh9EkdhEh70qYJWEREaqgJQkQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFElIBFRBJRAhYRSUQJWEQkESVgEZFE/h/c57DeUKuRMAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% Example with Frobenius norm + entropic regularization with gcg\n",
    "\n",
    "\n",
    "def f(G):\n",
    "    return 0.5 * np.sum(G**2)\n",
    "\n",
    "\n",
    "def df(G):\n",
    "    return G\n",
    "\n",
    "\n",
    "reg1 = 1e-3\n",
    "reg2 = 1e-1\n",
    "\n",
    "Gel2 = ot.optim.gcg(a, b, M, reg1, reg2, f, df, verbose=True)\n",
    "\n",
    "pl.figure(5, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, Gel2, 'OT entropic + matrix Frob. reg')\n",
    "pl.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
